﻿using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using System.Windows.Forms;
using System.Drawing;

namespace GeoFly
{
    public class RouteFlux
    {
        /// <summary>
        /// 栅格上的入流量
        /// </summary>
        public double InFlux;
        /// <summary>
        /// 存放一个栅格上的汇流输出量
        /// </summary>
        public double OutFlux;
        public bool bCal;
        public RouteFlux()
        {
            InFlux = 0.0;
            OutFlux = 0.0;
            bCal = false;
        }
    }
    //水文模拟总控制类
    public class HydroSimulate
    {
        public HydroSimulate()
        {
        }
        public HydroSimulate(DEMRiverNet river)
        {
            this.River = river;
        }
        public DEMRiverNet River;
        public static GridLayer GridDem = new GridLayer();
        public static GridLayerPara g_GridLayerPara = new GridLayerPara();
        /// <summary>
        /// 气象参数对象
        /// </summary>
        public static ClimatePara g_ClimatePara = new ClimatePara();
        public static MidGridResultOut g_GridResultOut = new MidGridResultOut();
        public static ModelRunPara g_ModelRunPara = new ModelRunPara();
        public static ModelInputPara g_ModelInputPara = new ModelInputPara();
        public static string path="C:\\ESSI_Data";
        /// <summary>
        /// 汇流次序表
        /// </summary>
        //public List<Cell> River.CalcuOrderList = new List<Cell>();


        /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
        +				径流过程模拟公用函数定义				+
        +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        +				控制：									+
        +					1、数据读入；						+
        +					2、内存初始化；						+
        +					3、结果输出；						+
        +++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
        /// <summary>
        /// 读入栅格图层参数（DEM,LULC,SOIL），该函数已修改
        /// </summary>
        /// <returns></returns>
        public void ReadInGridLayerData()
        {
            //读入DEM数据
            string filename = g_GridLayerPara.DemFileName;
            g_GridLayerPara.g_DemLayer = new GridLayer();
            g_GridLayerPara.g_DemLayer.ReadASC(filename);

            //读入土地利用数据
            filename = g_GridLayerPara.LulcFileName;
            g_GridLayerPara.g_VegLayer = new GridLayer();
            g_GridLayerPara.g_VegLayer.ReadASC(filename);

            //读入土壤类型数据
            filename = g_GridLayerPara.SoilFileName;
            g_GridLayerPara.g_SoilLayer = new GridLayer();
            g_GridLayerPara.g_SoilLayer.ReadASC(filename);

        }
        /// <summary>
        /// 读入气象数据参数
        /// </summary>
        /// <returns></returns>
        public bool ReadInClimateData()
        {
            MapExtent extent = g_GridLayerPara.g_DemLayer.Extent;
            double resolution = g_GridLayerPara.g_DemLayer.resolution;
            g_ClimatePara.pcpdata.SetExtentInfo(extent, resolution);
            g_ClimatePara.slrdata.SetExtentInfo(extent, resolution);
            g_ClimatePara.hmddata.SetExtentInfo(extent, resolution);
            g_ClimatePara.tmpmxdata.SetExtentInfo(extent, resolution);
            g_ClimatePara.tmpmndata.SetExtentInfo(extent, resolution);
            g_ClimatePara.tmpmeandata.SetExtentInfo(extent, resolution);
            g_ClimatePara.WindData.SetExtentInfo(extent, resolution);
            g_ClimatePara.petdata.SetExtentInfo(extent, resolution);

            //读入站点位置数据和值
            g_ClimatePara.pcpdata.ReadInStationPos(g_ClimatePara.PcpInfo.StationFileName);
            g_ClimatePara.pcpdata.ReadInStationValue(g_ClimatePara.PcpInfo.DataFileName);
            int count = g_ClimatePara.pcpdata.saDate.Count;
            DateTime startDate = g_ClimatePara.pcpdata.saDate[0];  //第一天
            DateTime endDate = g_ClimatePara.pcpdata.saDate[count - 1];//最后一天


            g_ClimatePara.StartDate = startDate;
            g_ClimatePara.EndDate = endDate;
            //如果蒸散发有现成的实际测量值，则读入实际测量值
            if (g_ModelRunPara.PETMethod == PET_METHOD.PET_REAL)
            {
                g_ClimatePara.petdata.ReadInStationPos(g_ClimatePara.PetInfo.StationFileName);
                g_ClimatePara.petdata.ReadInStationValue(g_ClimatePara.PetInfo.DataFileName);
            }

            //如果为暴雨
            if (g_ModelRunPara.RunoffSimuType == RunOffSimuType.STORM_RUNOFF_SIMULATION)
                return true;

            if (g_ClimatePara.SlrInfo.StationFileName != null && g_ClimatePara.SlrInfo.StationFileName != "")
            {
                g_ClimatePara.slrdata.ReadInStationPos(g_ClimatePara.SlrInfo.StationFileName);
                g_ClimatePara.slrdata.ReadInStationValue(g_ClimatePara.SlrInfo.DataFileName);
            }

            if (g_ClimatePara.HmdInfo.StationFileName != null && g_ClimatePara.HmdInfo.StationFileName != "")
            {
                g_ClimatePara.hmddata.ReadInStationPos(g_ClimatePara.HmdInfo.StationFileName);
                g_ClimatePara.hmddata.ReadInStationValue(g_ClimatePara.HmdInfo.DataFileName);
            }

            if (g_ClimatePara.WndInfo.StationFileName != null && g_ClimatePara.WndInfo.StationFileName!="")
            {
                g_ClimatePara.WindData.ReadInStationPos(g_ClimatePara.WndInfo.StationFileName);
                g_ClimatePara.WindData.ReadInStationValue(g_ClimatePara.WndInfo.DataFileName);
            }

            if (g_ClimatePara.TempmxInfo.StationFileName != null && g_ClimatePara.TempmxInfo.StationFileName != "")
            {
                g_ClimatePara.tmpmxdata.ReadInStationPos(g_ClimatePara.TempmxInfo.StationFileName);
                g_ClimatePara.tmpmxdata.ReadInStationValue(g_ClimatePara.TempmxInfo.DataFileName);
                g_ClimatePara.tmpmndata.ReadInStationPos(g_ClimatePara.TempmxInfo.StationFileName);
                g_ClimatePara.tmpmndata.ReadInStationValue(g_ClimatePara.TempmnInfo.DataFileName);

                g_ClimatePara.tmpmeandata.ReadInStationPos(g_ClimatePara.TempmxInfo.StationFileName);
                if (g_ClimatePara.TempmeanInfo.DataFileName != null && g_ClimatePara.TempmeanInfo.DataFileName !="")
                {
                    g_ClimatePara.tmpmeandata.ReadInStationValue(g_ClimatePara.TempmeanInfo.DataFileName);
                }
            }
            
            return true;
        }
        /// <summary>
        /// 读入马斯京根汇流参数
        /// </summary>
        /// <returns></returns>
        public bool ReadMuskingCoeff(string m_MuskCoeffFile)
        {
            if (File.Exists(m_MuskCoeffFile))
            {
                MessageBox.Show("马斯京根法（先演后合）汇流文件 [muskingum_coeff.txt] 不存在。");
                return false;
            }
            List<string> saValue = MatrixFuncs.FileRead(m_MuskCoeffFile);

            this.m_iNodeNum = Convert.ToInt32(saValue[0]);
            if (m_iNodeNum > 1)
            {
                m_pX = new double[m_iNodeNum];
                m_pK = new double[m_iNodeNum];

                int id = 0;
                for (int i = 2; i < saValue.Count; i++)
                {
                    string[] saOut = saValue[i].Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);

                    m_pX[id] = Convert.ToDouble(saOut[1]);
                    m_pK[id] = Convert.ToDouble(saOut[2]);
                    id++;
                }
            }
            return true;
        }

        public bool ReadWaterYearType(string WaterYearTypeFile)			//读入流域丰、平、枯水年份信息
        {
            List<string> saValue = MatrixFuncs.FileRead(WaterYearTypeFile);

            int num = saValue.Count;
            if (num > 0)
            {
                waterYearType = new SortedList<int, WaterYearType>();
                for (int i = 1; i < num; i++)
                {
                    string[] saOut = saValue[i].Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
                    int key = Convert.ToInt32(saOut[0]);
                    waterYearType[key] = (WaterYearType)Convert.ToInt32(saOut[1]);
                }
                return true;
            }
            else
                return false;
        }

        /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
        +			长时段降雨～径流过程模拟函数定义			+
        +++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
        /// <summary>
        /// 长时段降雨径流模拟主函数
        /// </summary>
        public void LongTermRunoffSimulate()
        {
            if (g_ModelRunPara.SurfRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
            {
                //读入栅格汇流最优次序参数文件
                this.ReadInRoutingPara();
            }

            this.ReadInRoutingLayerData();

            if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_COMBINE_FIRST ||
                g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_ROUTE_FIRST)
            {
                OpenFileDialog dlg = new OpenFileDialog();
                if (dlg.ShowDialog() == DialogResult.OK)
                {
                    ReadMuskingCoeff(dlg.FileName);
                }
            }

            int m_row = g_GridLayerPara.g_DemLayer.rowCount;
            int m_col = g_GridLayerPara.g_DemLayer.colCount;


            int totrec = 0; //= GetTimeInterval(true);//true表明按整日算
            this.MyOutletQ = new SortedList<DateTime, OutletQ>();


            if (m_iNodeNum > 1)
            {
                this.m_pNodeSurfQ = new SortedList<DateTime, double[]>();// new double[totrec, g_GridLayerPara.m_subNum];// GridOutletInit(totrec, g_GridLayerPara.m_subNum, 0);
                this.m_pNodeLatQ = new SortedList<DateTime, double[]>();// new double[totrec, g_GridLayerPara.m_subNum]; //GridOutletInit(totrec, g_GridLayerPara.m_subNum, 0);
                this.m_pNodeBaseQ = new SortedList<DateTime, double[]>();// new double[totrec, g_GridLayerPara.m_subNum]; //GridOutletInit(totrec, g_GridLayerPara.m_subNum, 0);
                this.m_pNodeOutQ = new SortedList<DateTime, double[]>();// new double[totrec, g_GridLayerPara.m_subNum]; //GridOutletInit(totrec, g_GridLayerPara.m_subNum, 0);
            }


            double dthet1334, dCp;
            double dsnowfactor = 1.0;
            double aetfactor = 0;

            int TotalYears = 0;
            int dn;

            dCp = 0.5;

            int StartYear, EndYear, startday, endday/*,mon,day*/;


            StartYear = g_ClimatePara.StartDate.Year;
            EndYear = g_ClimatePara.EndDate.Year;
            startday = DateAndTime.GetDnInYear(g_ClimatePara.StartDate.Year, g_ClimatePara.StartDate.Month, g_ClimatePara.StartDate.Day);
            endday = DateAndTime.GetDnInYear(g_ClimatePara.EndDate.Year, g_ClimatePara.EndDate.Month, g_ClimatePara.EndDate.Day);

            DateTime curorder = new DateTime(StartYear, g_ClimatePara.StartDate.Month, startday);

            TotalYears = EndYear - StartYear + 1;
            if (!ReadWaterYearType(this.m_WaterYrTypeFile))
            {
                waterYearType = new SortedList<int, WaterYearType>();
                for (int i = 0; i < TotalYears; i++)
                {
                    int key = StartYear + i;
                    waterYearType[key] = WaterYearType.WATER_LOW_YEAR;
                }
            }

            //水文过程年循环
            for (int year = StartYear; year <= EndYear; year++)
            {
                //按年份取得该年的模拟时段
                int ibeg = 0;
                int iend = 1;
                this.GetBEDateInYear(year, StartYear, EndYear, startday, endday, ibeg, iend);
                //从该年的模拟时段的开始日到结束日
                for (int j = ibeg; j <= iend; j++)
                {
                    dn = j;
                    PcpDataInterp(curorder);
                    MeteoDataInterp(curorder);


                    int iMonth = DateAndTime.GetMonthByDn(year, j);

                    double dintensity = 0;
                    double dhr = 24;
                    double dhrIntensity = 0;
                    double dPE = 0;
                    //对每个网格做下面的处理
                    for (int row = 0; row < m_row; row++)
                    {
                        for (int col = 0; col < m_col; col++)
                        {   //如果该网格已不需要计算，则跳过
                            if (!IfGridNeedCal(row, col))
                                continue;

                            dsnowfactor = 1;
                            //如果已读取了平均气温则执行

                            //如果温度低于雪温，则将降水加入雪量中
                            if (g_ClimatePara.tmpmeandata.m_gridLayer[row, col] < g_ModelInputPara.SnowTemperature)
                            {
                                g_GridResultOut.m_SnowWater[row, col] += g_ClimatePara.pcpdata.m_gridLayer[row, col];
                                //降雨量设为0
                                g_ClimatePara.pcpdata.m_gridLayer[row, col] = 0;

                                dsnowfactor = 0.15;
                            }
                            else	//温度高于雪温时
                            {
                                if (g_GridResultOut.m_SnowWater[row, col] > 0) //如果雪量大于0
                                {
                                    //求融雪量
                                    double smelt = DDFSnowMelt(g_ClimatePara.tmpmeandata.m_gridLayer[row, col],
                                                                g_ModelInputPara.SnowTemperature, g_ModelInputPara.DDF,
                                                                g_ModelInputPara.DailyMeanPcpTime);
                                    //如果实际雪量小于融雪量，就把融雪量设为实际雪量
                                    if (g_GridResultOut.m_SnowWater[row, col] < smelt)
                                    {
                                        smelt = g_GridResultOut.m_SnowWater[row, col];
                                        g_GridResultOut.m_SnowWater[row, col] = 0;//雪量设为0
                                    }
                                    else
                                        g_GridResultOut.m_SnowWater[row, col] -= smelt;
                                    //降水量中加入融雪量
                                    g_ClimatePara.pcpdata.m_gridLayer[row, col] += smelt;
                                    dsnowfactor = 0.3;
                                }
                            }

                            //每日平均降水时间(小时数)	
                            dhrIntensity = g_ModelInputPara.DailyMeanPcpTime;
                            //每小时降水量
                            dintensity = g_ClimatePara.pcpdata.m_gridLayer[row, col] / dhrIntensity;
                            //设置网格参数
                            HortonInfil hortonInfil = new HortonInfil();
                            hortonInfil.SetGridPara(g_GridLayerPara, row, col, g_GridResultOut.m_SoilProfileWater.Values[row,col]);
                            //计算Horton超渗产流量
                            
                            //土壤水下渗量
                            m_drateinf[row, col] = hortonInfil.HortonExcessRunoff();
                            //计算本栅格上的土壤水赤字
                            int SoilID = (int)g_GridLayerPara.g_SoilLayer[row, col];
                            dthet1334 = g_GridLayerPara.SoilTypes[SoilID].SoilWaterDeficitContent(g_GridResultOut.m_SoilProfileWater.Values[row,col]);

                            GridWaterBalance gridwb = new GridWaterBalance();
                            gridwb.SetGridPara(row, col, dintensity,m_drateinf[row,col], curorder, dhrIntensity);
                            int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
                            double dalb = GetVegAlbedo(VegID, iMonth, 1);//获取当前月的反照率

                            gridwb.CalcPET(g_ModelRunPara, g_ClimatePara, dalb);

                            if (g_GridLayerPara.g_StrahlerRivNet[row, col] != 0)
                            {
                                g_GridResultOut.m_GridSurfQ[row, col] = g_ClimatePara.pcpdata.m_gridLayer[row, col] - gridwb.m_dPET;
                                if (g_GridResultOut.m_GridSurfQ[row, col] < 0)
                                    g_GridResultOut.m_GridSurfQ[row, col] = 0;
                                if (g_GridResultOut.m_GridSurfQ[row, col] > 1e+10)
                                    MessageBox.Show("有错误");//AfxMessageBox("hello");
                                g_GridResultOut.m_GridLateralQ[row, col] = 0;
                                g_GridResultOut.m_GridBaseQ[row, col] = 0;
                                g_GridResultOut.m_GridTotalQ[row, col] = g_GridResultOut.m_GridSurfQ[row, col];
                                g_GridResultOut.m_AET[row, col] = gridwb.m_dPET;
                                g_GridResultOut.m_CI[row, col] = 0;
                                g_GridResultOut.m_CIDefict[row, col] = 0;
                                g_GridResultOut.m_NetPcp[row, col] = g_GridResultOut.m_GridSurfQ[row, col];

                                g_GridResultOut.m_GridWaterYieldType[row, col] = 0;
                                g_GridResultOut.m_SoilProfileWater[row, col] = 0;
                                g_GridResultOut.m_SoilAvgWater[row, col] = 0;

                            }
                            else
                            {
                                gridwb.CalcCI();
                                g_GridResultOut.m_CI[row, col] = gridwb.m_dCrownInterc;

                                if (g_GridResultOut.m_SoilProfileWater.Values[row,col] / g_GridLayerPara.SoilTypes[SoilID].SP_Fc > 0.8)
                                {
                                    if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
                                        aetfactor = 0.6;
                                    else
                                        aetfactor = 0.9;
                                }
                                else
                                {
                                    if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
                                        aetfactor = 0.4;
                                    else
                                        aetfactor = 0.6;
                                }

                                //*************对蒸散发处理的特殊代码段 -- 计算实际蒸散发**************//
                                if (g_ModelRunPara.PETMethod == PET_METHOD.PET_REAL)
                                {
                                    gridwb.m_dAET = g_ClimatePara.petdata.m_gridLayer[row, col] * aetfactor;
                                    g_GridResultOut.m_AET[row, col] = gridwb.m_dAET;
                                }
                                else
                                {
                                    gridwb.CalcAET(g_ModelRunPara, g_ClimatePara, dalb);
                                    if (gridwb.m_dAET > g_ClimatePara.petdata.m_gridLayer[row, col])		// FRZZ MODI 2005-8-18
                                        g_GridResultOut.m_AET[row, col] = g_ClimatePara.petdata.m_gridLayer[row, col] *
                                                                                        Math.Exp(-1 * g_ClimatePara.petdata.m_gridLayer[row, col] / gridwb.m_dAET);
                                    if (gridwb.m_dAET < 0)
                                        g_GridResultOut.m_AET[row, col] = g_ClimatePara.petdata.m_gridLayer[row, col] * 0.1;
                                }
                                //**********************对蒸散发处理的特殊代码段**********************//


                                gridwb.CalcNetRain();
                                g_GridResultOut.m_CIDefict[row, col] = gridwb.m_dCIDeficit;
                                g_GridResultOut.m_NetPcp[row, col] = gridwb.m_dNetRain;

                                g_GridResultOut.m_GridWaterYieldType[row, col] = gridwb.CalcRunoffElement(g_GridLayerPara, g_ModelInputPara,g_GridResultOut.m_SoilProfileWater[row, col]);

                                g_GridResultOut.m_SoilProfileWater[row, col] = gridwb.Cell_SP_Sw; //g_GridLayerPara.SoilTypes[SoilID].SP_Sw;
                                g_GridResultOut.m_SoilAvgWater[row, col] = g_GridLayerPara.SoilTypes[SoilID].SoilAvgWater(gridwb.Cell_SP_Sw);
                                g_GridResultOut.m_GridTotalQ[row, col] = gridwb.m_dTotalQ;
                                g_GridResultOut.m_GridSurfQ[row, col] = gridwb.m_dSurfQ;
                                g_GridResultOut.m_GridLateralQ[row, col] = gridwb.m_dLateralQ;
                                g_GridResultOut.m_GridBaseQ[row, col] = gridwb.m_dBaseQ;

                                if (g_GridResultOut.m_GridSurfQ[row, col] > 1e+10)
                                    MessageBox.Show("hello2");
                            }
                        }
                    }

                    if (m_iNodeNum == 1 || g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_PURE_LAG)
                    {
                        PureLagGridRouting(g_GridResultOut.m_GridSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curorder, dsnowfactor, waterYearType[year]);
                        PureLagGridRouting(g_GridResultOut.m_GridLateralQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curorder, dsnowfactor, waterYearType[year]);
                        PureLagGridRouting(g_GridResultOut.m_GridBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curorder, dsnowfactor, waterYearType[year]);
                        this.MyOutletQ[curorder].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);
                        this.MyOutletQ[curorder].CalcOutletQ();
                    }
                    else
                    {
                        PureLagGridRouting_Node(g_GridResultOut.m_GridSurfQ, m_pNodeSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curorder, totrec, dsnowfactor, waterYearType[year]);
                        PureLagGridRouting_Node(g_GridResultOut.m_GridLateralQ, m_pNodeLatQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curorder, totrec, dsnowfactor, waterYearType[year]);
                        PureLagGridRouting_Node(g_GridResultOut.m_GridBaseQ, m_pNodeBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curorder, totrec, dsnowfactor, waterYearType[year - StartYear]);

                        double dSurf, dLat, dBase;
                        dSurf = dLat = dBase = 0;
                        //for(int i=0;i<g_GridLayerPara.m_subNum;i++)
                        //{
                        //    dSurf += m_pNodeSurfQ[curorder,i] * g_ModelInputPara1.SurfQLinearFactor;
                        //    dLat += m_pNodeLatQ[curorder,i];
                        //    dBase += m_pNodeBaseQ[curorder,i];
                        //}
                        this.MyOutletQ[curorder].m_pOutletSurfQ = dSurf;
                        this.MyOutletQ[curorder].m_pOutletLatQ = dLat;
                        this.MyOutletQ[curorder].m_pOutletBaseQ = dBase;
                        this.MyOutletQ[curorder].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);

                        if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_COMBINE_FIRST)
                        {
                            for (int i = 0; i < g_GridLayerPara.m_subNum; i++)
                            {
                                m_pNodeOutQ[curorder][i] = m_pNodeSurfQ[curorder][i] * g_ModelInputPara.SurfQLinearFactor
                                    + m_pNodeLatQ[curorder][i] + m_pNodeBaseQ[curorder][i];
                            }
                            MuskingumRiverRouting(24, m_pNodeOutQ, m_pX, m_pK, g_GridLayerPara.m_subNum, curorder);
                            this.MyOutletQ[curorder].m_pOutletQ = pRiverRoute[g_GridLayerPara.m_subNum - 1].pRoute.OutFlux
                                                        + this.MyOutletQ[curorder].m_pOutletDeepBaseQ;
                        }
                        else if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_ROUTE_FIRST)
                        {
                            ;
                        }
                        else
                        {
                            this.MyOutletQ[curorder].CalcOutletQ();
                        }
                    }
                    TimeSpan DaySpan = new TimeSpan(1, 0, 0, 0);
                    curorder = curorder + DaySpan;
                    DateTime dateTime = new DateTime(year, iMonth, j);
                    MidGridResultOut(dateTime, true);
                }
            }
            RiverOutletQ(true);
            m_iTotalRec = totrec;

        }
        /// <summary>
        /// 计算栅格上的snowfactor
        /// </summary>
        /// <param name="row"></param>
        /// <param name="col"></param>
        /// <returns></returns>
        private double CalcSnowFactor(int row, int col)
        {
            double SnowFactor = 1;
            //如果温度低于雪温，则将降水加入雪量中
            if (g_ClimatePara.tmpmeandata.m_gridLayer[row, col] < g_ModelInputPara.SnowTemperature)
            {
                g_GridResultOut.m_SnowWater[row, col] += g_ClimatePara.pcpdata.m_gridLayer[row, col];
                //降雨量设为0
                g_ClimatePara.pcpdata.m_gridLayer[row, col] = 0;

                SnowFactor = 0.15;
            }
            else	//温度高于雪温时
            {
                if (g_GridResultOut.m_SnowWater[row, col] > 0) //如果雪量大于0
                {
                    //求融雪量
                    double smelt = DDFSnowMelt(g_ClimatePara.tmpmeandata.m_gridLayer[row, col],
                                                g_ModelInputPara.SnowTemperature, g_ModelInputPara.DDF,
                                                g_ModelInputPara.DailyMeanPcpTime);
                    //如果实际雪量小于融雪量，就把融雪量设为实际雪量
                    if (g_GridResultOut.m_SnowWater[row, col] < smelt)
                    {
                        smelt = g_GridResultOut.m_SnowWater[row, col];
                        g_GridResultOut.m_SnowWater[row, col] = 0;//雪量设为0
                    }
                    else
                        g_GridResultOut.m_SnowWater[row, col] -= smelt;
                    //降水量中加入融雪量
                    g_ClimatePara.pcpdata.m_gridLayer[row, col] += smelt;
                    SnowFactor = 0.3;
                }
            }
            return SnowFactor;
        }
        /// <summary>
        /// 长时段降雨径流模拟主函数(刘永和修改)
        /// </summary>
        public void LongTermRunoffSimulate_Liu()
        {
            if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_COMBINE_FIRST ||
                g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_ROUTE_FIRST)
            {
                OpenFileDialog dlg = new OpenFileDialog();
                if (dlg.ShowDialog() == DialogResult.OK)
                {
                    ReadMuskingCoeff(dlg.FileName);
                }
            }

            int rowCount = g_GridLayerPara.g_DemLayer.rowCount;
            int colCount = g_GridLayerPara.g_DemLayer.colCount;


            int totrec = 0; //= GetTimeInterval(true);//true表明按整日算
            this.MyOutletQ = new SortedList<DateTime, OutletQ>();


            if (m_iNodeNum > 1)
            {
                this.m_pNodeSurfQ = new SortedList<DateTime, double[]>();
                this.m_pNodeLatQ = new SortedList<DateTime, double[]>();
                this.m_pNodeBaseQ = new SortedList<DateTime, double[]>();
                this.m_pNodeOutQ = new SortedList<DateTime, double[]>();
            }



            int TotalYears = g_ClimatePara.StartDate.Year - g_ClimatePara.EndDate.Year + 1;

            if (!this.ReadWaterYearType(this.m_WaterYrTypeFile))
            {
                waterYearType = new SortedList<int, WaterYearType>();
                for (int year = g_ClimatePara.StartDate.Year; year <= g_ClimatePara.EndDate.Year; year++)
                {
                    waterYearType[year] = WaterYearType.WATER_LOW_YEAR;
                }
            }


            //水文过程年循环
            for (DateTime curDate = g_ClimatePara.StartDate; curDate <= g_ClimatePara.EndDate; curDate.Add(new TimeSpan(1, 0, 0, 0)))
            {

                PcpDataInterp(curDate);
                MeteoDataInterp(curDate);


                //int iMonth = curDate.Month;

                double dintensity = 0;
                double dhr = 24;
                double dhrIntensity = 0;
                double dPE = 0;
                double dsnowfactor = 1.0;
                //对每个网格做下面的处理
                for (int row = 0; row < rowCount; row++)
                {
                    for (int col = 0; col < colCount; col++)
                    {   //如果该网格已不需要计算，则跳过
                        if (!IfGridNeedCal(row, col))
                            continue;

                        dsnowfactor = this.CalcSnowFactor(row, col);

                        //每日平均降水时间(小时数)	
                        dhrIntensity = g_ModelInputPara.DailyMeanPcpTime;
                        //每小时降水量
                        dintensity = g_ClimatePara.pcpdata.m_gridLayer[row, col] / dhrIntensity;

                        HortonInfil hortonInfil = new HortonInfil();
                        //设置网格参数
                        hortonInfil.SetGridPara(g_GridLayerPara, row, col,g_GridResultOut.m_SoilProfileWater.Values[row,col]);
                        //土壤水下渗量
                        m_drateinf[row, col] = hortonInfil.HortonExcessRunoff();
                        //计算本栅格上的土壤水赤字
                        int SoilID = (int)g_GridLayerPara.g_SoilLayer[row, col];
                        double dthet = g_GridLayerPara.SoilTypes[SoilID].SoilWaterDeficitContent(g_GridResultOut.m_SoilProfileWater.Values[row,col]);


                        GridWaterBalance gridwb = new GridWaterBalance();
                        gridwb.SetGridPara(row, col, dintensity,m_drateinf.Values[row,col], curDate, dhrIntensity);
                        int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
                        double dalb = GetVegAlbedo(VegID, curDate.Month, 1);//获取当前月的反照率
                        double PET = gridwb.CalcPET(g_ModelRunPara, g_ClimatePara, dalb);

                        //如果当前栅格为河道
                        if (g_GridLayerPara.g_StrahlerRivNet[row, col] != 0)
                        {
                            g_GridResultOut.m_GridSurfQ[row, col] = g_ClimatePara.pcpdata.m_gridLayer[row, col] - gridwb.m_dPET;
                            if (g_GridResultOut.m_GridSurfQ[row, col] < 0)
                                g_GridResultOut.m_GridSurfQ[row, col] = 0;
                            if (g_GridResultOut.m_GridSurfQ[row, col] > 1e+10)
                                MessageBox.Show("有错误");
                            g_GridResultOut.m_GridLateralQ[row, col] = 0;
                            g_GridResultOut.m_GridBaseQ[row, col] = 0;
                            g_GridResultOut.m_GridTotalQ[row, col] = g_GridResultOut.m_GridSurfQ[row, col];
                            g_GridResultOut.m_AET[row, col] = PET;
                            g_GridResultOut.m_CI[row, col] = 0;
                            g_GridResultOut.m_CIDefict[row, col] = 0;
                            g_GridResultOut.m_NetPcp[row, col] = g_GridResultOut.m_GridSurfQ[row, col];

                            g_GridResultOut.m_GridWaterYieldType[row, col] = 0;
                            g_GridResultOut.m_SoilProfileWater[row, col] = 0;
                            g_GridResultOut.m_SoilAvgWater[row, col] = 0;

                        }
                        else //当前栅格不是河道
                        {
                            g_GridResultOut.m_CI[row, col] = gridwb.CalcCI();
                            double AetFactor = 0;
                            if (g_GridResultOut.m_SoilProfileWater.Values[row,col] / g_GridLayerPara.SoilTypes[SoilID].SP_Fc > 0.8)
                            {
                                if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
                                    AetFactor = 0.6;
                                else
                                    AetFactor = 0.9;
                            }
                            else
                            {
                                if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
                                    AetFactor = 0.4;
                                else
                                    AetFactor = 0.6;
                            }

                            //*************对蒸散发处理的特殊代码段 -- 计算实际蒸散发**************//
                            if (g_ModelRunPara.PETMethod == PET_METHOD.PET_REAL)
                            {
                                gridwb.m_dAET = g_ClimatePara.petdata.m_gridLayer[row, col] * AetFactor;
                                g_GridResultOut.m_AET[row, col] = gridwb.m_dAET;
                            }
                            else
                            {
                                gridwb.CalcAET(g_ModelRunPara, g_ClimatePara, dalb);
                                if (gridwb.m_dAET > g_ClimatePara.petdata.m_gridLayer[row, col])		// FRZZ MODI 2005-8-18
                                    g_GridResultOut.m_AET[row, col] = g_ClimatePara.petdata.m_gridLayer[row, col] * Math.Exp(-1 * g_ClimatePara.petdata.m_gridLayer[row, col] / gridwb.m_dAET);
                                if (gridwb.m_dAET < 0)
                                    g_GridResultOut.m_AET[row, col] = g_ClimatePara.petdata.m_gridLayer[row, col] * 0.1;
                            }
                            //**********************对蒸散发处理的特殊代码段**********************//


                            gridwb.CalcNetRain();
                            g_GridResultOut.m_CIDefict[row, col] = gridwb.m_dCIDeficit;
                            g_GridResultOut.m_NetPcp[row, col] = gridwb.m_dNetRain;

                            g_GridResultOut.m_GridWaterYieldType[row, col] = gridwb.CalcRunoffElement(g_GridLayerPara, g_ModelInputPara,g_GridResultOut.m_SoilProfileWater[row, col]);
                            g_GridResultOut.m_SoilProfileWater[row, col] = gridwb.Cell_SP_Sw; //g_GridLayerPara.SoilTypes[SoilID].SP_Sw;
                            g_GridResultOut.m_SoilAvgWater[row, col] = g_GridLayerPara.SoilTypes[SoilID].SoilAvgWater(gridwb.Cell_SP_Sw);
                            g_GridResultOut.m_GridTotalQ[row, col] = gridwb.m_dTotalQ;
                            g_GridResultOut.m_GridSurfQ[row, col] = gridwb.m_dSurfQ;
                            g_GridResultOut.m_GridLateralQ[row, col] = gridwb.m_dLateralQ;
                            g_GridResultOut.m_GridBaseQ[row, col] = gridwb.m_dBaseQ;

                            if (g_GridResultOut.m_GridSurfQ[row, col] > 1e+10)
                                MessageBox.Show("hello2");
                        }
                    }
                }

                if (m_iNodeNum == 1 || g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_PURE_LAG)
                {
                    PureLagGridRouting(g_GridResultOut.m_GridSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curDate, dsnowfactor, waterYearType[curDate.Year]);
                    PureLagGridRouting(g_GridResultOut.m_GridLateralQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curDate, dsnowfactor, waterYearType[curDate.Year]);
                    PureLagGridRouting(g_GridResultOut.m_GridBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curDate, dsnowfactor, waterYearType[curDate.Year]);
                    int dn = DateAndTime.GetDnInYear(curDate.Year, curDate.Month, curDate.Day);
                    this.MyOutletQ[curDate].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);
                    this.MyOutletQ[curDate].m_pOutletQ = this.MyOutletQ[curDate].m_pOutletSurfQ * g_ModelInputPara.SurfQLinearFactor
                        + this.MyOutletQ[curDate].m_pOutletLatQ + this.MyOutletQ[curDate].m_pOutletBaseQ + this.MyOutletQ[curDate].m_pOutletDeepBaseQ;
                }
                else
                {
                    PureLagGridRouting_Node(g_GridResultOut.m_GridSurfQ, m_pNodeSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curDate, totrec, dsnowfactor, waterYearType[curDate.Year]);
                    PureLagGridRouting_Node(g_GridResultOut.m_GridLateralQ, m_pNodeLatQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curDate, totrec, dsnowfactor, waterYearType[curDate.Year]);
                    PureLagGridRouting_Node(g_GridResultOut.m_GridBaseQ, m_pNodeBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curDate, totrec, dsnowfactor, waterYearType[curDate.Year]);

                    double dSurf, dLat, dBase;
                    dSurf = dLat = dBase = 0;
                    for (int i = 0; i < g_GridLayerPara.m_subNum; i++)
                    {
                        //dSurf += m_pNodeSurfQ[curorder, i] * g_ModelInputPara1.SurfQLinearFactor;
                        //dLat += m_pNodeLatQ[curorder, i];
                        //dBase += m_pNodeBaseQ[curorder, i];
                    }
                    this.MyOutletQ[curDate].m_pOutletSurfQ = dSurf;
                    this.MyOutletQ[curDate].m_pOutletLatQ = dLat;
                    this.MyOutletQ[curDate].m_pOutletBaseQ = dBase;
                    int dn = DateAndTime.GetDnInYear(curDate.Year, curDate.Month, curDate.Day);
                    this.MyOutletQ[curDate].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);

                    if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_COMBINE_FIRST)
                    {
                        for (int i = 0; i < g_GridLayerPara.m_subNum; i++)
                        {
                            m_pNodeOutQ[curDate][i] = m_pNodeSurfQ[curDate][i] * g_ModelInputPara.SurfQLinearFactor
                                + m_pNodeLatQ[curDate][i] + m_pNodeBaseQ[curDate][i];
                        }
                        MuskingumRiverRouting(24, m_pNodeOutQ, m_pX, m_pK, g_GridLayerPara.m_subNum, curDate);
                        this.MyOutletQ[curDate].m_pOutletQ = pRiverRoute[g_GridLayerPara.m_subNum - 1].pRoute.OutFlux
                                                    + this.MyOutletQ[curDate].m_pOutletDeepBaseQ;
                    }
                    else if (g_ModelRunPara.RiverRouteMethod == RiverRouteMethod.ROUTE_MUSKINGUM_ROUTE_FIRST)
                    {
                        ;
                    }
                    else
                    {
                        this.MyOutletQ[curDate].m_pOutletQ = this.MyOutletQ[curDate].m_pOutletSurfQ +
                                                    this.MyOutletQ[curDate].m_pOutletLatQ +
                                                    this.MyOutletQ[curDate].m_pOutletBaseQ +
                                                    this.MyOutletQ[curDate].m_pOutletDeepBaseQ;
                    }
                }

                MidGridResultOut(curDate, true);
            }

            RiverOutletQ(true);
            m_iTotalRec = totrec;

        }

        /// <summary>
        /// 通过标准日期类型转化为字符串(Julian历格式)
        /// the result maybe like 1998001,1998002,...1998365,
        /// </summary>
        /// <param name="date">日期</param>
        /// <returns>格式字符串</returns>
        public string SetLongDateByDn(DateTime date)
        {
            int dn = DateAndTime.GetDnInYear(date.Year, date.Month, date.Day);
            return dn.ToString();
        }





        /// <summary>
        /// 根据年份和Julian历日期，返回日期字符串
        /// </summary>
        /// <param name="year">年份（4位）</param>
        /// <param name="dn">dn -- Julian历日期</param>
        /// <returns>日期字符串（格式：1998-1-1）</returns>
        public DateTime DnToDate(int year, int dn)
        {
            int date = 0;
            int[] DnNum = new int[12];
            for (int i = 1; i <= 12; i++)
                DnNum[i - 1] = DateAndTime.GetEachMonthDn(year, i);
            DateTime datetime = new DateTime();
            if (dn <= DnNum[0])
            {
                datetime = new DateTime(year, 1, dn);
            }
            else
            {
                for (int i = 1; i < 12; i++) //other month
                {
                    if (dn <= DnNum[i] && dn > DnNum[i - 1])
                    {
                        date = dn - DnNum[i - 1];
                        datetime = new DateTime(year, i + 1, date);
                    }
                }
            }
            return datetime;
        }



        /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
        +				暴雨～径流过程模拟函数定义				+
        +++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
        
        /// <summary>
        /// 单独计算一个栅格的水量平衡
        /// </summary>
        /// <param name="row">行号</param>
        /// <param name="col">列号</param>
        /// <param name="mon">月份</param>
        /// <param name="dhr">当前时段的小时数</param>
        /// <param name="curorder">当前执行序号</param>
        /// <param name="dn">年内的日序号</param>
        public void CalcCell_Horton(int row, int col, double dhr, DateTime curorder, int dn)
        {
            //如果该网格点不需要计算，则跳过
            if (this.IfGridNeedCal(row, col) == false)
                return;
            int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
            double dalb = this.GetVegAlbedo(VegID, curorder.Month, 1);
            //如果栅格点上有河流级别，执行下面
            if (River.m_pStrahlerOrd[row, col] != 0)
            {
                //降水量减去蒸散发，取得净降水量
                if (g_ClimatePara.pcpdata.m_gridLayer[row, col] < 0 || g_ClimatePara.petdata.m_gridLayer[row, col] < 0)
                    throw new Exception("降雨或蒸发量数据不能为负值");
                g_GridResultOut.m_GridSurfQ[row, col] = g_ClimatePara.pcpdata.m_gridLayer[row, col] - g_ClimatePara.petdata.m_gridLayer[row, col];
                //如果得到的差为负值，则地表径流设为零
                if (g_GridResultOut.m_GridSurfQ[row, col] < 0)
                {
                    g_GridResultOut.m_GridSurfQ[row, col] = 0;
                }
                //栅格壤中流和栅格地下径流都设为0(这是因为河道内可以说无壤中流和潜流)
                g_GridResultOut.m_GridLateralQ[row, col] = 0;
                g_GridResultOut.m_GridBaseQ[row, col] = 0;
            }
            else//如果该栅格无河流级别
            {
                //计算降雨强度(以每小时计算(mm/h))
                double dintensity = g_ClimatePara.pcpdata.m_gridLayer[row, col] / dhr;
                //计算Horton超渗产流
                HortonInfil hortonInfil= new HortonInfil();
                hortonInfil.SetGridPara(g_GridLayerPara, row, col, g_GridResultOut.m_SoilProfileWater[row,col]);
                //计算本栅格土壤水实际下渗量
                m_drateinf[row, col] = hortonInfil.HortonExcessRunoff();

                
                
                //下面采用水平衡类进行计算
                GridWaterBalance gridwb = new GridWaterBalance();
                //设置栅格上用于计算水平衡的参数
                gridwb.SetGridPara(row, col, dintensity,m_drateinf.Values[row,col], curorder, dhr);

                double aetfactor = 1.0;
                int SoilID = (int)g_GridLayerPara.g_SoilLayer[row, col];
                if (g_GridResultOut.m_SoilProfileWater[row,col] / g_GridLayerPara.SoilTypes[SoilID].SP_Fc > 0.8)
                {
                    if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
                        aetfactor = 0.75;
                    else
                        aetfactor = 0.95;
                }
                else
                {
                    if (g_ClimatePara.pcpdata.m_gridLayer[row, col] > 0)
                        aetfactor = 0.45;
                    else
                        aetfactor = 0.65;
                }

                //*************对蒸散发处理的特殊代码段 -- 计算实际蒸散发**************//                            
                gridwb.m_dAET = g_ClimatePara.petdata.m_gridLayer[row, col] * aetfactor;
                g_GridResultOut.m_AET[row, col] = gridwb.m_dAET;

                //计算净雨量
                gridwb.CalcNetRain();
                g_GridResultOut.m_NetPcp[row, col] = gridwb.m_dNetRain;
                g_GridResultOut.m_GridWaterYieldType[row, col] = gridwb.CalcRunoffElement(g_GridLayerPara, g_ModelInputPara,g_GridResultOut.m_SoilProfileWater[row, col]);

                g_GridResultOut.m_SoilProfileWater[row, col] = gridwb.Cell_SP_Sw; 
                g_GridResultOut.m_SoilAvgWater[row, col] = g_GridLayerPara.SoilTypes[SoilID].SoilAvgWater(gridwb.Cell_SP_Sw);
                g_GridResultOut.m_GridTotalQ[row, col] = gridwb.m_dTotalQ;
                g_GridResultOut.m_GridSurfQ[row, col] = gridwb.m_dSurfQ;
                g_GridResultOut.m_GridLateralQ[row, col] = gridwb.m_dLateralQ;
                g_GridResultOut.m_GridBaseQ[row, col] = gridwb.m_dBaseQ;
                //Console.WriteLine(gridwb.m_dLateralQ);
            }
        }
        public void CalcCell_GreenAmpt(int row, int col, double dhr, DateTime curorder, int dn)
        {
            int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
            double dalb = this.GetVegAlbedo(VegID, curorder.Month, 1);

            string key = g_ClimatePara.petdata.saStationName[0];

            //**********************对蒸散发处理的特殊代码段**********************//
            //如果当前栅格的strahler分级为河流
            if (River.m_pStrahlerOrd[row, col] != 0)
            {
                g_GridResultOut.m_GridSurfQ[row, col] = g_ClimatePara.pcpdata.m_gridLayer[row, col] - g_ClimatePara.petdata.m_gridLayer[row, col];
                if (g_GridResultOut.m_GridSurfQ[row, col] < 0.0)
                    g_GridResultOut.m_GridSurfQ[row, col] = 0.0;
                g_GridResultOut.m_GridLateralQ[row, col] = 0.0;
                g_GridResultOut.m_GridBaseQ[row, col] = 0.0;
            }
            else //strahler分级不是河流
            {
                double  dintensity = g_ClimatePara.pcpdata.m_gridLayer[row, col] / dhr;//降雨强度
                double dCp = g_GridLayerPara.VegTypes[VegID].CoverDeg[curorder.Month - 1];
                GreenAmptInfil GreenAmpt = new GreenAmptInfil(g_GridLayerPara);
                GreenAmpt.SetGridPara(row, col, 1.0, m_drateinf[row, col], m_dcumr[row, col], m_drintns[row, col], m_dcuminf[row, col], m_dexcum[row, col], g_ClimatePara.pcpdata.m_gridLayer[row, col], dCp);

                GreenAmpt.GreenAmptExcessRunoff(g_GridResultOut.m_SoilProfileWater.Values[row, col]);
                m_drateinf[row, col] = GreenAmpt.m_dPreRateinf;
                m_dcumr[row, col] = GreenAmpt.m_dPrecumr;
                m_drintns[row, col] = GreenAmpt.m_dPrerintns;
                m_dcuminf[row, col] = GreenAmpt.m_dPrecuminf;
                m_dexcum[row, col] = GreenAmpt.m_dPreexcum;
                int SoilID = (int)g_GridLayerPara.g_SoilLayer[row, col];                

                GridWaterBalance gridwb = new GridWaterBalance();
                gridwb.SetGridPara(row, col, dintensity, m_drateinf.Values[row, col], curorder, dhr);

                
                //**********************对蒸散发处理的特殊代码段**********************//
                if (g_ClimatePara.PetInfo.StationFileName == "" && g_ClimatePara.PetInfo.DataFileName != "")
                {
                    key = g_ClimatePara.petdata.saStationName[0];
                    gridwb.m_dAET = g_ClimatePara.petdata.pMeteoData[key][curorder];
                    g_GridResultOut.m_AET[row, col] = gridwb.m_dAET;
                }
                //**********************对蒸散发处理的特殊代码段**********************//
                
                gridwb.CalcNetRain();
                g_GridResultOut.m_NetPcp[row, col] = gridwb.m_dNetRain;

                g_GridResultOut.m_GridWaterYieldType[row, col] = gridwb.CalcRunoffElement(g_GridLayerPara, g_ModelInputPara, g_GridResultOut.m_SoilProfileWater[row, col]);
                g_GridResultOut.m_SoilProfileWater[row, col] = gridwb.Cell_SP_Sw;
                g_GridResultOut.m_SoilAvgWater[row, col] = g_GridLayerPara.SoilTypes[SoilID].SoilAvgWater(gridwb.Cell_SP_Sw);
                g_GridResultOut.m_GridTotalQ[row, col] = gridwb.m_dTotalQ;
                g_GridResultOut.m_GridSurfQ[row, col] = gridwb.m_dSurfQ;
                g_GridResultOut.m_GridLateralQ[row, col] = gridwb.m_dLateralQ;
                g_GridResultOut.m_GridBaseQ[row, col] = gridwb.m_dBaseQ;
            }
        }
        /// <summary>
        /// 进行汇流计算
        /// </summary>
        /// <param name="curorder">当前的执行时段序号</param>
        /// <param name="dhr">当前时段的小时数</param>
        /// <param name="totrec">总时段数</param>
        /// <param name="dn">年内的日序</param>
        public void CalcRiverRouting_Horton(DateTime curorder, double dhr, int dn)
        {
            //下面进行坡面漫流汇流
            if (g_ModelRunPara.SurfRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
            {
                this.MuskCungeGridRouting(dhr, g_GridResultOut.m_GridSurfQ, g_ModelInputPara.SMCTimeWeight, g_ModelInputPara.SMCGridRTravelTime);
                for (int k = 0; k < RoutePara.Length; k++)
                {
                    int row = RoutePara[k].pRow;
                    int col = RoutePara[k].pCol;
                    g_GridResultOut.m_GridRoutingQ[row, col] = pSurfQ[k].pRoute.OutFlux;
                }
                //把路径上最后一个像素的出流量作为总出流量
                if (this.MyOutletQ.ContainsKey(curorder) == false)
                    this.MyOutletQ[curorder] = new OutletQ();
                this.MyOutletQ[curorder].m_pOutletSurfQ = pSurfQ[RoutePara.Length - 1].pRoute.OutFlux;
            }
            else
            {
                this.PureLagGridRouting(g_GridResultOut.m_GridSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curorder);
            }
            //下面进行壤中流汇流
            if (g_ModelRunPara.LatRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
            {
                MuskCungeGridRouting(dhr, g_GridResultOut.m_GridLateralQ, g_ModelInputPara.LMCTimeWeight, g_ModelInputPara.LMCGridRTravelTime);

                for (int k = 0; k < RoutePara.Length; k++)
                {
                    g_GridResultOut.m_GridRoutingQ.Values[RoutePara[k].pRow, RoutePara[k].pCol] += pLatQ[k].pRoute.OutFlux;
                }
                if (this.MyOutletQ.ContainsKey(curorder) == false)
                    this.MyOutletQ[curorder] = new OutletQ();
                this.MyOutletQ[curorder].m_pOutletLatQ = pLatQ[RoutePara.Length - 1].pRoute.OutFlux;
            }
            else
            {
                PureLagGridRouting(g_GridResultOut.m_GridLateralQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curorder);
            }
            //下面进行地下径流的汇流
            if (g_ModelRunPara.BaseRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
            {
                this.MuskCungeGridRouting(dhr, g_GridResultOut.m_GridBaseQ, g_ModelInputPara.BMCTimeWeight, g_ModelInputPara.BMCGridRTravelTime);

                for (int k = 0; k < RoutePara.Length; k++)
                {
                    g_GridResultOut.m_GridRoutingQ.Values[RoutePara[k].pRow, RoutePara[k].pCol] += pBaseQ[k].pRoute.OutFlux;
                }
                if (this.MyOutletQ.ContainsKey(curorder) == false)
                    this.MyOutletQ[curorder] = new OutletQ();
                this.MyOutletQ[curorder].m_pOutletBaseQ = pBaseQ[RoutePara.Length - 1].pRoute.OutFlux;
            }
            else
            {
                PureLagGridRouting(g_GridResultOut.m_GridBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curorder);
            }
            //下面计算深层地下径流的汇流
            this.MyOutletQ[curorder].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);
            this.MyOutletQ[curorder].m_pOutletQ = this.MyOutletQ[curorder].m_pOutletSurfQ * g_ModelInputPara.SurfQLinearFactor
                                        + this.MyOutletQ[curorder].m_pOutletLatQ + this.MyOutletQ[curorder].m_pOutletBaseQ + this.MyOutletQ[curorder].m_pOutletDeepBaseQ;
        }
        
        /// <summary>
        /// 格林－安普特下渗曲线计算暴雨过程主函数
        /// </summary>
        /// <returns></returns>
        public void StormRunoffSim_GreenAmpt_Liu()
        {
            int rowCount = g_GridLayerPara.g_DemLayer.rowCount;
            int colCount = g_GridLayerPara.g_DemLayer.colCount;
            this.CellOutSurfFlux = new double[rowCount, colCount];
            this.CellOutLatFlux = new double[rowCount, colCount];
            this.CellOutBaseFlux = new double[rowCount, colCount];
            GridLayerInit();
            GridLayerInit_GreenAmpt();
            
            this.MyOutletQ = new SortedList<DateTime, OutletQ>();
            int n = River.CalcuOrderList.Count;
            CurrentSurfQ = new RouteFlux[n];
            PreSurfQ = new RouteFlux[n];
            for (int i = 0; i < n; i++)
            {
                CurrentSurfQ[i] = new RouteFlux();
                PreSurfQ[i] = new RouteFlux();
            }
            StreamWriter sw = new StreamWriter(HydroSimulate.path + "\\Output\\OutletQ_Liu.txt");
            sw.WriteLine("Order\tTotalQ\tSurfQ\tLatQ\tBaseQ\tDeepBaseQ");

            //对各降雨时段进行循环
            int count = g_ClimatePara.pcpdata.saDate.Count;
            for (int i = 0; i < count; i++)
            {
                //功能：给定年、月、日，计算该日的Julian历日
                
                DateTime curorder = g_ClimatePara.pcpdata.saDate[i];
                //对i时段的降雨量进行内插
                MeteoDataInterp(curorder);	//气象数据内插
                //间隔的小时数
                double dhr = 1.0;
                int dn = DateAndTime.GetDnInYear(curorder.Year, curorder.Month, curorder.Day);

                for (int row = 0; row < rowCount; row++)
                {
                    for (int col = 0; col < colCount; col++)
                    {
                        if (!IfGridNeedCal(row, col))
                            continue;
                        this.CalcCell_GreenAmpt(row, col, 1.0, curorder, dn);                        
                    }
                }

                this.CalcRiverRouting_Liu(curorder, 1, dn);

                sw.Write(curorder + "\t");
                sw.Write(this.MyOutletQ[curorder].m_pOutletQ.ToString("0.0000") + "\t");
                sw.Write(this.MyOutletQ[curorder].m_pOutletSurfQ.ToString("0.0000") + "\t");
                sw.Write(this.MyOutletQ[curorder].m_pOutletLatQ.ToString("0.0000") + "\t");
                sw.Write(this.MyOutletQ[curorder].m_pOutletBaseQ.ToString("0.0000") + "\t");
                sw.Write(this.MyOutletQ[curorder].m_pOutletDeepBaseQ.ToString("0.0000") + "\t");
                sw.WriteLine();
                sw.Flush();
                MidGridResultOut(g_ClimatePara.pcpdata.saDate[i], false);
            }
            //输出断面流量结果
            this.RiverOutletQ(true);
            // true;
        }
        /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
        +					流域汇流模拟函数定义				+
        +++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
        /// <summary>
        /// 读入栅格汇流最优次序参数文件
        /// </summary>
        /// <returns></returns>
        public void ReadInRoutingPara()
        {
            m_OutRow = m_OutCol = 0;

            string strFileStub = g_GridLayerPara.DemFileName.Substring(0, g_GridLayerPara.DemFileName.Length - 4);
            string strFileName = strFileStub + "_gud.txt";
            if (!File.Exists(strFileName))
                throw new Exception("栅格汇流文件" + strFileName + "不存在");

            List<string> saIn = MatrixFuncs.FileRead(strFileName);

            int LineCount = saIn.Count; //行数
            RoutePara = new MuskCungeRoutePara[LineCount];
            pRouteQ = new double[LineCount];
            pSurfQ = new MuskRouteFlux[LineCount];
            pLatQ = new MuskRouteFlux[LineCount];
            pBaseQ = new MuskRouteFlux[LineCount];


            ProgressBar bar = new ProgressBar();
            bar.Text = "正在加载全栅格汇流参数文件";
            bar.Show();
            for (int row = 0; row < LineCount; row++)
            {
                string strLine = saIn[row];
                string[] saLine = strLine.Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
                RoutePara[row] = new MuskCungeRoutePara();
                RoutePara[row].pGridNum = Convert.ToInt32(saLine[0]);
                RoutePara[row].pGridRouteOrd = Convert.ToInt32(saLine[1]);
                for (int i = 0; i < 8; i++)
                {
                    RoutePara[row].pInGrid[i] = Convert.ToInt32(saLine[i + 2]);
                }
                RoutePara[row].pOutGrid = Convert.ToInt32(saLine[10]);
                RoutePara[row].pGridSlope = Convert.ToDouble(saLine[11]);
                RoutePara[row].pGridRLength = Convert.ToDouble(saLine[12]);
                RoutePara[row].pGridRiverOrd = Convert.ToInt32(saLine[13]);
                RoutePara[row].pRow = Convert.ToInt32(saLine[14]);
                RoutePara[row].pCol = Convert.ToInt32(saLine[15]);

                pSurfQ[row] = new MuskRouteFlux();
                pSurfQ[row].pRoute = new RouteFlux();
                pSurfQ[row].pPreRoute = new RouteFlux();
                pSurfQ[row].pPreRoute.bCal = true;

                pLatQ[row] = new MuskRouteFlux();
                pLatQ[row].pRoute = new RouteFlux();
                pLatQ[row].pPreRoute = new RouteFlux();
                pLatQ[row].pPreRoute.bCal = true;

                pBaseQ[row] = new MuskRouteFlux();
                pBaseQ[row].pRoute = new RouteFlux();
                pBaseQ[row].pPreRoute = new RouteFlux();
                pBaseQ[row].pPreRoute.bCal = true;

                bar.progressBar1.Value = (int)(row * 1.0 / LineCount * 100);
            }
            bar.Close();
            m_OutRow = RoutePara[LineCount - 1].pRow;
            m_OutCol = RoutePara[LineCount - 1].pCol;
        }

        /// <summary>
        /// 读入栅格汇流相关参数，包括全流域边界、子流域边界、河道等级矩阵、地表径流汇流参数矩阵
        /// 壤中径流汇流参数矩阵、地下径流汇流参数矩阵、地下径流参数矩阵
        /// </summary>
        /// <returns></returns>
        public void ReadInRoutingLayerData()
        {
            string strFileStub = g_GridLayerPara.DemFileName.Substring(0, g_GridLayerPara.DemFileName.Length - 4);

            string strFileName = strFileStub + "_wby.asc";
            g_GridLayerPara.g_BasinBoundary = new GridLayer();
            g_GridLayerPara.g_BasinBoundary.ReadASC(strFileName);


            strFileName = strFileStub + "_sws.asc";
            g_GridLayerPara.g_SubWaterShed = new GridLayer();
            g_GridLayerPara.g_SubWaterShed.ReadASC(strFileName);
            g_GridLayerPara.g_SubWaterShed.CalMaxMinValue();
            g_GridLayerPara.m_subNum = (int)g_GridLayerPara.g_SubWaterShed.MaxValue;


            strFileName = strFileStub + "_sor.asc";
            g_GridLayerPara.g_StrahlerRivNet = new GridLayer();
            g_GridLayerPara.g_StrahlerRivNet.ReadASC(strFileName);
            g_GridLayerPara.g_StrahlerRivNet.CalMaxMinValue();
            g_GridLayerPara.m_MaxStrahlerOrd = (int)g_GridLayerPara.g_StrahlerRivNet.MaxValue;


            if (g_ModelRunPara.SurfRouteMethod == SurfRouteMehtod.ROUTE_PURE_LAG)
            {
                strFileName = strFileStub + "_gst.asc";
                g_GridLayerPara.g_RouteSurfQTime.ReadASC(strFileName);
            }

            if (g_ModelRunPara.LatRouteMethod == SurfRouteMehtod.ROUTE_PURE_LAG)
            {
                strFileName = strFileStub + "_glt.asc";
                g_GridLayerPara.g_RouteLatQTime = new GridLayer();
                g_GridLayerPara.g_RouteLatQTime.ReadASC(strFileName);
            }

            if (g_ModelRunPara.BaseRouteMethod == SurfRouteMehtod.ROUTE_PURE_LAG)
            {
                strFileName = strFileStub + "_gbt.asc";
                g_GridLayerPara.g_RouteBaseQTime = new GridLayer();
                g_GridLayerPara.g_RouteBaseQTime.ReadASC(strFileName);
            }
        }


        /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
        +				径流过程模拟私有函数定义				+
        +++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

        public void GridLayerInit()
        {
            GridLayer DemLayer=g_GridLayerPara.g_DemLayer;

            g_GridResultOut.m_GridTotalQ = DemLayer.AttributesCopy();
            g_GridResultOut.m_GridSurfQ = DemLayer.AttributesCopy();
            g_GridResultOut.m_GridLateralQ = DemLayer.AttributesCopy();
            g_GridResultOut.m_GridBaseQ = DemLayer.AttributesCopy();            
            g_GridResultOut.m_GridWaterYieldType = DemLayer.AttributesCopy();
            
            g_GridResultOut.m_SoilAvgWater = DemLayer.AttributesCopy();
            g_GridResultOut.m_GridRoutingQ = DemLayer.AttributesCopy();
            g_GridResultOut.m_NetPcp = DemLayer.AttributesCopy();
            g_GridResultOut.m_AET = DemLayer.AttributesCopy();
            g_GridResultOut.m_AI = DemLayer.AttributesCopy();
            g_GridResultOut.m_CI = DemLayer.AttributesCopy();
            g_GridResultOut.m_CIDefict = DemLayer.AttributesCopy();
            g_GridResultOut.m_NetPcp = DemLayer.AttributesCopy();
            g_GridResultOut.m_GridWaterYieldType = DemLayer.AttributesCopy();
            m_drateinf = DemLayer.AttributesCopy(0);
            g_GridResultOut.m_SoilProfileWater = DemLayer.AttributesCopy();
            for (int row = 0; row < DemLayer.rowCount; row++)
            {
                for (int col = 0; col < DemLayer.colCount; col++)
                {
                    if (River.m_pGridRouteLength[row, col] > 0)
                    {
                        int SoilID =(int)g_GridLayerPara.g_SoilLayer[row, col];
                        SoilInfo soilInfo = g_GridLayerPara.SoilTypes[SoilID];
                        g_GridResultOut.m_SoilProfileWater[row, col] = soilInfo.Initial_SP_Sw;
                    }                        
                }
            }
            
        }



        /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
        +														+
        +					气象数据空间插值		    		+
        +														+
        +++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
        /// <summary>
        /// 降雨量空间插值（有时段性问题，所以和其它气象数据空间插值分开）
        /// </summary>
        /// <param name="ord"></param>
        public void PcpDataInterp(DateTime timeOrd)
        {
            g_ClimatePara.pcpdata.SpatialGridOut(timeOrd);
        }

        public void MeteoDataInterp(DateTime timeOrd)
        {
            
            GridLayer pcp= g_ClimatePara.pcpdata.SpatialGridOut(timeOrd);
            //pcp.WriteASC(path + "\\Test\\pcp.txt");
            //g_ClimatePara.tmpmxdata.SpatialGridOut(timeOrd, 2, false);
            //g_ClimatePara.tmpmndata.SpatialGridOut(timeOrd, 2, false);
            //g_ClimatePara.slrdata.SpatialGridOut(timeOrd, 2, false);
            //g_ClimatePara.hmddata.SpatialGridOut(timeOrd, 2, false);
            //g_ClimatePara.WindData.SpatialGridOut(timeOrd, 2, false);
            g_ClimatePara.petdata.SpatialGridOut(timeOrd);
            //g_ClimatePara.tmpmeandata.SpatialGridOut(timeOrd, 2, false);
            //else
            //{
            //    if (g_ClimatePara.bReadTmpMX == true && g_ClimatePara.bReadTmpMN == true)
            //    {
            //        g_ClimatePara.tmpmeandata.m_gridLayer.extent = g_ClimatePara.tmpmndata.m_gridLayer.extent.Copy();
            //        g_ClimatePara.tmpmeandata.m_gridLayer.resolution = g_ClimatePara.tmpmndata.m_gridLayer.resolution;

            //        g_ClimatePara.tmpmeandata.m_gridLayer.initGridStruct();
            //        for (int j = 0; j < g_ClimatePara.tmpmeandata.m_gridLayer.rowCount; j++)
            //        {
            //            for (int k = 0; k < g_ClimatePara.tmpmeandata.m_gridLayer.colCount; k++)
            //            {
            //                double dvalue;
            //                dvalue = (g_ClimatePara.tmpmxdata.m_gridLayer.Values[j, k] + g_ClimatePara.tmpmndata.m_gridLayer.Values[j, k]) / 2.0;
            //                g_ClimatePara.tmpmeandata.m_gridLayer.Values[j, k] = dvalue;
            //            }
            //        }
            //    }
            //}
        }
        /// <summary>
        /// 网格中的中间结果输出
        /// </summary>
        /// <param name="yr"></param>
        /// <param name="dn"></param>
        /// <param name="ord"></param>
        /// <param name="bDate"></param>
        public void MidGridResultOut(DateTime dateTime, bool bDate)
        {

            DateTime strCurDate = dateTime;
            string TimeString = strCurDate.Year + "_" + strCurDate.Month + "_" + strCurDate.Day + "_" + strCurDate.Hour;
            string strtmp;



            if (HydroSimulate.g_GridResultOut.iPcp)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\pcp_" + TimeString + ".asc";
                HydroSimulate.g_ClimatePara.pcpdata.m_gridLayer.WriteASC(strtmp);
            }
            if (g_GridResultOut.iTempMax)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\TempMax_" + TimeString + ".asc";
                HydroSimulate.g_ClimatePara.tmpmxdata.m_gridLayer.WriteASC(strtmp);
            }
            if (g_GridResultOut.iTempMin)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\TempMin_" + TimeString + ".asc";
                HydroSimulate.g_ClimatePara.tmpmndata.m_gridLayer.WriteASC(strtmp);
            }
            if (g_GridResultOut.iTempMean)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\TempMean_" + TimeString + ".asc";
                HydroSimulate.g_ClimatePara.tmpmeandata.m_gridLayer.WriteASC(strtmp);
            }
            if (g_GridResultOut.iSlr)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\Solar_" + TimeString + ".asc";
                HydroSimulate.g_ClimatePara.slrdata.m_gridLayer.WriteASC(strtmp);
            }
            if (g_GridResultOut.iWnd)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\Wnd_" + TimeString + ".asc";
                HydroSimulate.g_ClimatePara.WindData.m_gridLayer.WriteASC(strtmp);
            }
            if (g_GridResultOut.iHmd)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\Hmd_" + TimeString + ".asc";
                HydroSimulate.g_ClimatePara.hmddata.m_gridLayer.WriteASC(strtmp);
            }
            if (g_GridResultOut.iPET)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\PET_" + TimeString + ".asc";
                HydroSimulate.g_ClimatePara.petdata.m_gridLayer.WriteASC(strtmp);
            }
            if (g_GridResultOut.iAET)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\AET_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_AET.WriteASC(strtmp);
            }
            if (g_GridResultOut.iCI)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\CI_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_CI.WriteASC(strtmp);

            }
            if (g_GridResultOut.iSnowWater)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\Snow_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_SnowWater.WriteASC(strtmp);

            }
            if (g_GridResultOut.iAI)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\AI_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_AI.WriteASC(strtmp);
            }
            if (g_GridResultOut.iRouteOut)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\GridRouteQ_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_GridRoutingQ.WriteASC(strtmp);
            }
            if (g_GridResultOut.iSurfQ)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\SurfQ_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_GridSurfQ.WriteASC(strtmp);

            }
            if (g_GridResultOut.iWaterYieldType)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\WYT_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_GridWaterYieldType.WriteASC(strtmp);
            }
            if (g_GridResultOut.iLatQ)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\LatQ_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_GridLateralQ.WriteASC(strtmp);
            }
            if (g_GridResultOut.iBaseQ)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\BaseQ_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_GridBaseQ.WriteASC(strtmp);
            }
            if (g_GridResultOut.iInfilRate)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\InfilRate_" + TimeString + ".asc";
                HydroSimulate.m_drateinf.WriteASC(strtmp);
            }
            if (g_GridResultOut.iProfileSoilWater)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\SoilProfileWater_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_SoilProfileWater.WriteASC(strtmp);

            }
            if (g_GridResultOut.iAvgSoilWater)
            {
                strtmp = HydroSimulate.g_GridResultOut.OutputPath + "\\SoilAvgWater_" + TimeString + ".asc";
                HydroSimulate.g_GridResultOut.m_SoilAvgWater.WriteASC(strtmp);

            }

        }

        /// <summary>
        /// 输出所有流量值
        /// </summary>
        /// <param name="recnum">时段</param>
        /// <param name="bDate"></param>
        public void RiverOutletQ(bool bDate)
        {            
            if(!Directory.Exists(HydroSimulate.path+"\\Output"))
                Directory.CreateDirectory(HydroSimulate.path+"\\Output");
            StreamWriter sw = new StreamWriter(HydroSimulate.path + "\\Output\\OutQ.txt");
            sw.WriteLine("Order\tTotalQ\tSurfQ\tLatQ\tBaseQ\tDeepBaseQ");
            if (bDate)
            {
                foreach (DateTime timeKey in this.MyOutletQ.Keys)
                {
                    sw.Write(timeKey + "\t");
                    sw.Write(this.MyOutletQ[timeKey].m_pOutletQ.ToString("0.0000") + "\t");
                    sw.Write(this.MyOutletQ[timeKey].m_pOutletSurfQ.ToString("0.0000") + "\t");
                    sw.Write(this.MyOutletQ[timeKey].m_pOutletLatQ.ToString("0.0000") + "\t");
                    sw.Write(this.MyOutletQ[timeKey].m_pOutletBaseQ.ToString("0.0000") + "\t");
                    sw.Write(this.MyOutletQ[timeKey].m_pOutletDeepBaseQ.ToString("0.0000") + "\t");
                    sw.WriteLine();
                }
            }
            sw.Close();
        }
        /// <summary>
        /// 判断当前网格是否需要计算
        /// </summary>
        /// <param name="row"></param>
        /// <param name="col"></param>
        /// <returns>true-需要计算，false-不需要计算</returns>
        public bool IfGridNeedCal(int row, int col)
        {

            //如果该网格点是流域边界，则不需要计算
            //if (g_GridLayerPara.g_BasinBoundary[row, col] != 1)
            //    return false;
            if (this.River.m_pGridRouteLength[row, col] == 0)
                return false;
            //如果网格点上的植被数据为空
            if (g_GridLayerPara.g_VegLayer[row, col] == g_GridLayerPara.g_VegLayer.nodata)
                return false;
            //如果网格点上的土壤数据为空
            if (g_GridLayerPara.g_SoilLayer[row, col] == g_GridLayerPara.g_SoilLayer.nodata)
                return false;
            //如果网格点上的DEM数据为空
            if (g_GridLayerPara.g_DemLayer[row, col] == g_GridLayerPara.g_DemLayer.nodata)
                return false;
            
            //不符合以上条件的则需要计算
            return true;
        }

        /// <summary>
        /// 度－日因子模型计算日融雪量
        /// </summary>
        /// <param name="dtav"></param>
        /// <param name="dtThresh"></param>
        /// <param name="ddf"></param>
        /// <param name="dtlen"></param>
        /// <returns></returns>
        public double DDFSnowMelt(double dtav, double dtThresh, double ddf, double dtlen)
        {
            return ddf * (dtav - dtThresh) * dtlen;
        }

        public double DeepBaseQSim(int dn, double dmin)
        {
            double dret = 0;
            if (dmin <= 0)
                return dret;
            else
            {
                dret = (dmin + 3 * dmin) / 2 + (3 * dmin - dmin) / 2 * Math.Sin(2 * Math.PI * (dn - 82) / 365);
            }
            return dret;
        }

        //public void GridLayerInit_LongTerm()
        //{
        //    m_drateinf.initGridStruct();
        //    m_drintns.initGridStruct();
        //    m_dcumr.initGridStruct();
        //    m_dcuminf.initGridStruct();
        //    m_dexcum.initGridStruct();
        //    g_GridResultOut.m_AI.initGridStruct();
        //    g_GridResultOut.m_CI.initGridStruct();
        //    g_GridResultOut.m_SnowWater.initGridStruct();
        //    g_GridResultOut.m_CIDefict.initGridStruct();
        //    g_GridResultOut.m_AET.initGridStruct();
        //    g_ClimatePara.petdata.m_gridLayer.initGridStruct();
        //    g_GridResultOut.m_NetPcp.initGridStruct();
        //}
        
        /// <summary>
        /// 按年份取得该年的模拟时段
        /// </summary>
        /// <param name="yr">指定年份</param>
        /// <param name="startyr">起始年</param>
        /// <param name="endyr">结束年</param>
        /// <param name="startday">起始天</param>
        /// <param name="endday">结束天</param>
        /// <param name="ibeg">存储起始日</param>
        /// <param name="iend">存储结束日</param>
        public void GetBEDateInYear(int yr, int startyr, int endyr, int startday, int endday, int ibeg, int iend)
        {
            bool bLeap = DateAndTime.CheckLeapYear(yr);
            if (bLeap)
                iend = 366;
            else
                iend = 365;

            if (yr == startyr)
                ibeg = startday;
            else
                ibeg = 1;

            if (yr == endyr)
                iend = endday;
        }
        

        
       
        
        public void GridLayerInit_GreenAmpt()
        {
            

            m_drintns = g_GridLayerPara.g_DemLayer.AttributesCopy(0);

            m_dcumr = g_GridLayerPara.g_DemLayer.AttributesCopy(0);

            m_dcuminf = g_GridLayerPara.g_DemLayer.AttributesCopy(0);

            m_dexcum = g_GridLayerPara.g_DemLayer.AttributesCopy(0);
            g_GridResultOut.m_AI = g_GridLayerPara.g_DemLayer.AttributesCopy();

            g_GridResultOut.m_CI = g_GridLayerPara.g_DemLayer.AttributesCopy();

            g_GridResultOut.m_CIDefict = g_GridLayerPara.g_DemLayer.AttributesCopy();
            g_GridResultOut.m_AET = g_GridLayerPara.g_DemLayer.AttributesCopy();

            g_ClimatePara.petdata.m_gridLayer = g_GridLayerPara.g_DemLayer.AttributesCopy();
        }

        
        

        /// <summary>
        /// 马斯京根-康吉法汇流计算(先合后演)
        /// </summary>
        /// <param name="dT">时间的离散值</param>
        /// <param name="GridQ"></param>
        /// <param name="pRoute"></param>
        /// <param name="pPreRoute"></param>
        /// <param name="x"></param>
        /// <param name="k"></param>
        /// <returns></returns>
        public void MuskCungeGridRouting(double dT, GridLayer GridQ, double x, double k)
        {
            ProgressBar bar = new ProgressBar();
            bar.Text = "正在计算汇流";
            bar.Show();
            //按照汇流次序对逐个栅格进行汇流
            for (int i = 0; i < RoutePara.Length; i++)
            {
                //先统计当前栅格(即第i个)的总入流量
                double SumQ = 0;
                for (int j = 0; j < 8; j++)
                {
                    if (RoutePara[i].pInGrid[j] == 0) //第i个栅格在j方向上无入流，则绕过
                        continue;
                    if (pSurfQ[RoutePara[i].pInGrid[j] - 1].pRoute.bCal == true)
                        SumQ += pSurfQ[RoutePara[i].pInGrid[j] - 1].pRoute.OutFlux;
                    else
                        throw new Exception("routing grid uncalculated");

                }
                pSurfQ[i].pRoute.InFlux = SumQ;

                //再使用Cunge法计算出流量
                pSurfQ[i].pRoute.OutFlux = GridRoutingOut(i, RoutePara[i].pGridSlope, RoutePara[i].pGridRLength, dT,
                                                        pSurfQ[i].pPreRoute.InFlux, pSurfQ[i].pRoute.InFlux, pSurfQ[i].pPreRoute.OutFlux,
                                                        GridQ.Values[RoutePara[i].pRow, RoutePara[i].pCol] * Math.Pow(GridQ.resolution, 2) / (1000 * dT * 3600),//栅格面入流
                                                        x, k);
                pSurfQ[i].pRoute.bCal = true;
                bar.progressBar1.Value = (int)(i * 100.0 / RoutePara.Length);

            }
            bar.Close();
            //标识新的入流和出流关系
            for (int i = 0; i < RoutePara.Length; i++)
            {
                pSurfQ[i].pPreRoute.InFlux = pSurfQ[i].pRoute.InFlux;
                pSurfQ[i].pPreRoute.OutFlux = pSurfQ[i].pRoute.OutFlux;
                pSurfQ[i].pRoute.bCal = false;
            }

        }
        /// <summary>
        /// 计算序号为ord的栅格的汇流结果，采用Cunge法
        /// </summary>
        /// <param name="ord">栅格的序号</param>
        /// <param name="gridslp">网格坡度</param>
        /// <param name="deltaX">河长的离散值</param>
        /// <param name="deltaT">时间的离散值</param>
        /// <param name="Q11">差分值</param>
        /// <param name="Q12">差分值</param>
        /// <param name="Q21">差分值</param>
        /// <param name="dGridQ">当前时段的侧向流量</param>
        /// <param name="x">河段</param>
        /// <param name="k">时段</param>
        /// <returns>流出流量</returns>
        public double GridRoutingOut(int ord, double gridslp, double deltaX, double deltaT, double Q11, double Q12, double Q21, double dGridQ, double x, double k)
        {
            double dVflow = 0;
            RiverType iRiverType;
            double dB = 0;//用于存放河道或山的宽度
            double dUnitQ = 0;
            double Q22 = 0;
            //如果是河道，计算河道流速
            if (RoutePara[ord].pGridRiverOrd > 0)
            {
                iRiverType = g_ModelInputPara.RiverProType;
                dB = g_ModelInputPara.RiverProWidth;
                dUnitQ = dGridQ / dB;//当前时段的侧向单宽入流
                double slope = RoutePara[ord].pGridSlope;//取出坡度
                int row = RoutePara[ord].pRow;
                int col = RoutePara[ord].pCol;
                int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
                double dManning = g_GridLayerPara.VegTypes[VegID].dManning_n;
                dVflow = 0.489 * Math.Pow(dUnitQ, 0.25) * Math.Pow(slope, 0.375) / Math.Pow(dManning, 0.75);
            }
            else //如果是坡面，计算坡面流速
            {
                iRiverType = g_ModelInputPara.HillProType;
                dB = g_ModelInputPara.HillProWidth;
                if (dB == 0)
                    dB = g_GridLayerPara.g_DemLayer.resolution;
                //计算当前时段的侧向单宽入流
                dUnitQ = dGridQ / dB;
                //计算坡面流速
                int row = RoutePara[ord].pRow;
                int col = RoutePara[ord].pCol;
                int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
                double dManning = g_GridLayerPara.VegTypes[VegID].dManning_n;
                dVflow = Math.Pow(dUnitQ, 0.4) * Math.Pow(RoutePara[ord].pGridSlope, 0.3) / Math.Pow(dManning, 0.6);
            }
            //准备Cunge法所需要的参数
            MuskCunge.SetRoutingPara(dVflow, gridslp, deltaX, deltaT, iRiverType, Q11, Q12, Q21, dUnitQ, dB, x, k);
            //计算Cunge法的洪水波流量
            Q22 = MuskCunge.RoutingOutQ();
            return Q22;
        }
        /// <summary>
        /// 计算序号为ord的栅格的汇流结果，采用Cunge法
        /// </summary>
        /// <param name="ord">栅格的序号</param>
        /// <param name="gridslp">网格坡度</param>
        /// <param name="deltaX">河长的离散值</param>
        /// <param name="deltaT">时间的离散值</param>
        /// <param name="Q11">差分值</param>
        /// <param name="Q12">差分值</param>
        /// <param name="Q21">差分值</param>
        /// <param name="dGridQ">当前时段的侧向流量</param>
        /// <param name="TimeWeight">河段</param>
        /// <param name="TravelTime">时段</param>
        /// <returns>流出流量</returns>
        public double GridRoutingOut_Liu(int ord, double deltaX, double deltaT, double Q11, double Q12, double Q21, double dGridQ, double TimeWeight, double TravelTime)
        {
            double dVflow = 0;
            RiverType iRiverType;
            double dB = 0;//用于存放河道或山的宽度
            double dUnitQ = 0;

            //如果是河道，计算河道流速
            int row = River.CalcuOrderList[ord].row;
            int col = River.CalcuOrderList[ord].col;
            double slope = this.River.m_pSlope[row, col];//取出坡度
            double dManning = 0.01;
            int VegID = (int)g_GridLayerPara.g_VegLayer[row, col];
            if (VegID!=-9999 && g_GridLayerPara.VegTypes[VegID] != null)
                dManning = g_GridLayerPara.VegTypes[VegID].dManning_n;
            if (River.m_pStrahlerOrd[row, col] > 0)
            {
                iRiverType = g_ModelInputPara.RiverProType;
                dB = g_ModelInputPara.RiverProWidth;
                dUnitQ = dGridQ / dB;//当前时段的侧向单宽入流

                dVflow = 0.489 * Math.Pow(dUnitQ, 0.25) * Math.Pow(slope, 0.375) / Math.Pow(dManning, 0.75);
            }
            else //如果是坡面，计算坡面流速
            {
                iRiverType = g_ModelInputPara.HillProType;
                dB = g_ModelInputPara.HillProWidth;
                if (dB == 0)
                    dB = g_GridLayerPara.g_DemLayer.resolution;
                //计算当前时段的侧向单宽入流
                dUnitQ = dGridQ / dB;
                //计算坡面流速
                dVflow = Math.Pow(dUnitQ, 0.4) * Math.Pow(slope, 0.3) / Math.Pow(dManning, 0.6);
            }
            //准备Cunge法所需要的参数
            double x = TimeWeight;
            double k = TravelTime/60.0; //把分钟转化为小时
            double denominator = k * (1 - x) + 0.5 * deltaT;
            double C1 = (x * k + 0.5 * deltaT) / denominator;
            double C2 = (0.5 * deltaT - x * k) / denominator;
            double C3 = (k * (1 - x) - 0.5 * deltaT) / denominator;
            double C4 = deltaT * deltaX / denominator;
            //计算Cunge法的有限差分值
            double Q22= C1 * Q11 + C2 * Q12 + C3 * Q21 + C4 * dUnitQ;
            return Q22;
        }

        /*	马斯京根法(按子流域划分河段)  */

        public bool MuskingumRiverRouting(double dTLen, SortedList<DateTime, double[]> pNodeQ, double[] pX, double[] pK, int NodeNum, DateTime curOrder)
        {
            bool bret = true;
            for (int i = 0; i < NodeNum; i++)
            {
                double QOut = 0;
                if (i == 0)
                    QOut = pNodeQ[curOrder][i];
                else
                    QOut = pNodeQ[curOrder][i] + pRiverRoute[i - 1].pRoute.OutFlux;

                pRiverRoute[i].pRoute.InFlux = QOut;

                pRiverRoute[i].pRoute.OutFlux = RiverRoutingOut(dTLen, pRiverRoute[i].pPreRoute.InFlux, pRiverRoute[i].pRoute.InFlux,
                pRiverRoute[i].pPreRoute.OutFlux, pX[i], pK[i]);
                pRiverRoute[i].pRoute.bCal = true;

            }
            for (int i = 0; i < NodeNum; i++)
            {
                pRiverRoute[i].pPreRoute.InFlux = pRiverRoute[i].pRoute.InFlux;
                pRiverRoute[i].pPreRoute.OutFlux = pRiverRoute[i].pRoute.OutFlux;
                pRiverRoute[i].pRoute.bCal = false;
            }
            return bret;

        }


        public double RiverRoutingOut(double deltaT, double dInFlux1, double dInFlux2, double dOutFlux1, double x, double k)
        {
            double dOutFlux2 = 0;

            Muskingum.SetRoutingPara(deltaT, dInFlux1, dInFlux2, dOutFlux1, x, k);

            dOutFlux2 = Muskingum.RoutingOutQ();

            return dOutFlux2;
        }



        /// <summary>
        /// 滞时演算法(单流域)
        /// </summary>
        /// <param name="GridQ"></param>
        /// <param name="pRoute"></param>
        /// <param name="dTlen"></param>
        /// <param name="QType"></param>
        /// <param name="curorder"></param>
        /// <param name="totorder"></param>
        /// <returns></returns>
        public void PureLagGridRouting(GridLayer GridQ, double dTlen, RunOffSourceType QType, DateTime curorder)
        {
            this.PureLagGridRouting(GridQ, dTlen, QType, curorder, 1.0, WaterYearType.WATER_LOW_YEAR);

        }
        /// <summary>
        /// 滞时演算法(单流域)
        /// </summary>
        /// <param name="GridQ">待处理的数据</param>
        /// <param name="pRoute">路径像素集合</param>
        /// <param name="dTlen">时段按小时计算的长度</param>
        /// <param name="QType">径流的类型，如是壤中流，还是基流</param>
        /// <param name="curorder">当前时次</param>
        /// <param name="totorder"></param>
        /// <param name="snowfactor"></param>
        /// <param name="WaterYrType"></param>
        /// <returns></returns>
        public void PureLagGridRouting(GridLayer GridQ, double dhr, RunOffSourceType QType, DateTime curorder, double snowfactor, WaterYearType WaterYrType)
        {
            DateTime LagOrder;

            double dGridOut = 0;
            double dSurfQLoss = 0;


            if (WaterYrType == WaterYearType.WATER_HIGH_YEAR)
            {
                dSurfQLoss = g_ModelInputPara.HighWaterSurfQLoss;

            }
            else if (WaterYrType == WaterYearType.WATER_MID_YEAR)
            {
                dSurfQLoss = g_ModelInputPara.MidWaterSurfQLoss;

            }
            else if (WaterYrType == WaterYearType.WATER_LOW_YEAR)
            {
                dSurfQLoss = g_ModelInputPara.LowWaterSurfQLoss;
            }
            int m_row = g_GridLayerPara.g_DemLayer.rowCount;
            int m_col = g_GridLayerPara.g_DemLayer.colCount;
            for (int row = 0; row < m_row; row++)
            {
                for (int col = 0; col < m_col; col++)
                {
                    if (IfGridNeedCal(row, col) == false)
                        continue;
                    if (QType == RunOffSourceType.RUNOFF_ELEMENT_SURFQ)
                    {
                        int LagTime = (int)(g_GridLayerPara.g_RouteSurfQTime[row, col] / dhr); //滞后时间
                        TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
                        LagOrder = curorder.Add(LagSpan);//当前时段加上滞时
                        DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
                        if (LagOrder < endTime)
                        {
                            dGridOut = snowfactor * GridQ[row, col] / (dSurfQLoss * Math.Pow(g_GridLayerPara.g_RouteSurfQTime[row, col], 1));
                            if (this.MyOutletQ.ContainsKey(LagOrder))
                                this.MyOutletQ[LagOrder].m_pOutletSurfQ += dGridOut;
                        }
                    }
                    else if (QType == RunOffSourceType.RUNOFF_ELEMENT_LATERALQ)
                    {
                        int LagTime = (int)(g_GridLayerPara.g_RouteLatQTime[row, col] / dhr);
                        TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
                        LagOrder = curorder.Add(LagSpan);
                        DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
                        if (LagOrder < endTime)
                        {
                            dGridOut = snowfactor * GridQ[row, col] / (g_ModelInputPara.LatQLoss * Math.Pow(g_GridLayerPara.g_RouteLatQTime[row, col], 1));
                            if (this.MyOutletQ.ContainsKey(LagOrder))
                                this.MyOutletQ[LagOrder].m_pOutletLatQ += dGridOut;
                        }
                    }
                    else if (QType == RunOffSourceType.RUNOFF_ELEMENT_BASEQ)
                    {
                        int LagTime = (int)(g_GridLayerPara.g_RouteBaseQTime[row, col] / dhr);
                        TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
                        LagOrder = curorder.Add(LagSpan);
                        DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
                        if (LagOrder < endTime)
                        {
                            dGridOut = snowfactor * GridQ[row, col] / (g_ModelInputPara.BaseQLoss * Math.Pow(g_GridLayerPara.g_RouteBaseQTime[row, col], 1));
                            if (this.MyOutletQ.ContainsKey(LagOrder))
                                this.MyOutletQ[LagOrder].m_pOutletBaseQ += dGridOut;
                        }
                    }
                }

            }
        }
        /// <summary>
        /// 滞时演算法(单流域)
        /// </summary>
        /// <param name="GridQ">待处理的数据</param>
        /// <param name="pRoute">路径像素集合</param>
        /// <param name="dTlen">时段按小时计算的长度</param>
        /// <param name="QType">径流的类型，如是壤中流，还是基流</param>
        /// <param name="curorder">当前时次</param>
        /// <param name="totorder"></param>
        /// <param name="snowfactor"></param>
        /// <param name="WaterYrType"></param>
        /// <returns></returns>
        public void PureLagGridRouting_Liu(GridLayer GridQ, double dhr, RunOffSourceType QType, DateTime curorder, double snowfactor, WaterYearType WaterYrType)
        {
            DateTime LagOrder;

            double dGridOut = 0;
            double dSurfQLoss = 0;


            if (WaterYrType == WaterYearType.WATER_HIGH_YEAR)
            {
                dSurfQLoss = g_ModelInputPara.HighWaterSurfQLoss;

            }
            else if (WaterYrType == WaterYearType.WATER_MID_YEAR)
            {
                dSurfQLoss = g_ModelInputPara.MidWaterSurfQLoss;

            }
            else if (WaterYrType == WaterYearType.WATER_LOW_YEAR)
            {
                dSurfQLoss = g_ModelInputPara.LowWaterSurfQLoss;
            }
            int m_row = g_GridLayerPara.g_DemLayer.rowCount;
            int m_col = g_GridLayerPara.g_DemLayer.colCount;
            for (int row = 0; row < m_row; row++)
            {
                for (int col = 0; col < m_col; col++)
                {
                    if (IfGridNeedCal(row, col) == false)
                        continue;
                    if (QType == RunOffSourceType.RUNOFF_ELEMENT_SURFQ)
                    {
                        int LagTime = (int)(River.m_pSurfQRouteTime[row, col] / dhr); //滞后时间
                        TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
                        LagOrder = curorder.Add(LagSpan);//当前时段加上滞时
                        DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
                        if (LagOrder < endTime)
                        {
                            if (River.m_pSurfQRouteTime[row, col] <= 0)
                                continue;
                            dGridOut = snowfactor * GridQ[row, col] / (dSurfQLoss * Math.Pow(River.m_pSurfQRouteTime[row, col], 1));
                            if (this.MyOutletQ.ContainsKey(LagOrder))
                                this.MyOutletQ[LagOrder].m_pOutletSurfQ += dGridOut;
                            else
                            {
                                this.MyOutletQ[LagOrder] = new OutletQ();
                                this.MyOutletQ[LagOrder].m_pOutletSurfQ = dGridOut;
                            }
                        }
                    }
                    else if (QType == RunOffSourceType.RUNOFF_ELEMENT_LATERALQ)
                    {
                        int LagTime = (int)(River.m_pLatQRouteTime[row, col] / dhr);
                        TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
                        LagOrder = curorder.Add(LagSpan);
                        DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
                        if (LagOrder < endTime)
                        {
                            //无汇流
                            if (River.m_pLatQRouteTime[row, col] <= 0)
                                continue;
                            dGridOut = snowfactor * GridQ[row, col] / (g_ModelInputPara.LatQLoss * Math.Pow(River.m_pLatQRouteTime[row, col], 1));
                            //Console.Write(River.m_pBaseQRouteTime[row, col]+" ");
                            //Console.Write(g_ModelInputPara.LatQLoss+" ");
                            //Console.WriteLine(dGridOut);
                            if (this.MyOutletQ.ContainsKey(LagOrder))
                                this.MyOutletQ[LagOrder].m_pOutletLatQ += dGridOut;
                            else
                            {
                                this.MyOutletQ[LagOrder] = new OutletQ();
                                this.MyOutletQ[LagOrder].m_pOutletLatQ = dGridOut;
                            }
                        }
                    }
                    else if (QType == RunOffSourceType.RUNOFF_ELEMENT_BASEQ)
                    {
                        int LagTime = (int)(River.m_pBaseQRouteTime[row, col] / dhr);
                        TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
                        LagOrder = curorder.Add(LagSpan);
                        DateTime endTime = g_ClimatePara.pcpdata.saDate[g_ClimatePara.pcpdata.saDate.Count - 1];
                        if (LagOrder < endTime)
                        {
                            //无汇流
                            if (River.m_pBaseQRouteTime[row, col] <= 0)
                                continue;
                            dGridOut = snowfactor * GridQ[row, col] / (g_ModelInputPara.BaseQLoss * Math.Pow(River.m_pBaseQRouteTime[row, col], 1));
                            if (this.MyOutletQ.ContainsKey(LagOrder))
                                this.MyOutletQ[LagOrder].m_pOutletBaseQ += dGridOut;

                            else
                            {
                                this.MyOutletQ[LagOrder] = new OutletQ();
                                this.MyOutletQ[LagOrder].m_pOutletBaseQ = dGridOut;
                            }
                        }
                    }
                }

            }
        }

        /*	滞时演算法(分子流域)	*/
        public bool PureLagGridRouting_Node(GridLayer GridQ, SortedList<DateTime, double[]> pRoute, double dTlen,
                                RunOffSourceType QType,
                                DateTime curorder,
                                int totorder,
                                double snowfactor,
                                WaterYearType WaterYrType)
        {
            DateTime LagOrder;
            int LagTime = 0;
            double dGridOut = 0;
            int subID = 0;

            double dSurfQLoss = 0;

            switch (WaterYrType)
            {
                case WaterYearType.WATER_HIGH_YEAR:
                    {
                        dSurfQLoss = g_ModelInputPara.HighWaterSurfQLoss;
                        break;
                    }
                case WaterYearType.WATER_MID_YEAR:
                    {
                        dSurfQLoss = g_ModelInputPara.MidWaterSurfQLoss;
                        break;
                    }
                case WaterYearType.WATER_LOW_YEAR:
                    {
                        dSurfQLoss = g_ModelInputPara.LowWaterSurfQLoss;
                        break;
                    }
            }
            int m_row = g_GridLayerPara.g_DemLayer.rowCount;
            int m_col = g_GridLayerPara.g_DemLayer.colCount;
            for (int i = 0; i < m_row; i++)
            {
                for (int j = 0; j < m_col; j++)
                {
                    if (g_GridLayerPara.g_BasinBoundary.Values[i, j] == 1.0)
                    {
                        subID = (int)g_GridLayerPara.g_SubWaterShed.Values[i, j] - 1;
                        switch (QType)
                        {
                            case RunOffSourceType.RUNOFF_ELEMENT_SURFQ:
                                {
                                    LagTime = (int)(g_GridLayerPara.g_RouteSurfQTime.Values[i, j] / dTlen);
                                    TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
                                    LagOrder = curorder.Add(LagSpan);
                                    //if (LagOrder < totorder)
                                    //{
                                    if (g_GridLayerPara.g_RouteSurfQTime.Values[i, j] <= 0)
                                        g_GridLayerPara.g_RouteSurfQTime.Values[i, j] = 0.1;
                                    dGridOut = snowfactor * GridQ.Values[i, j] / (dSurfQLoss * Math.Pow(g_GridLayerPara.g_RouteSurfQTime.Values[i, j], 1));
                                    pRoute[LagOrder][subID] += dGridOut;
                                    //}
                                }
                                break;
                            case RunOffSourceType.RUNOFF_ELEMENT_LATERALQ:
                                {
                                    LagTime = (int)(g_GridLayerPara.g_RouteLatQTime.Values[i, j] / dTlen);
                                    TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
                                    LagOrder = curorder.Add(LagSpan);
                                    //if (LagOrder < totorder)
                                    //{
                                    if (g_GridLayerPara.g_RouteLatQTime.Values[i, j] <= 0)
                                        g_GridLayerPara.g_RouteLatQTime.Values[i, j] = 0.1;
                                    dGridOut = snowfactor * GridQ.Values[i, j] / (g_ModelInputPara.LatQLoss * Math.Pow(g_GridLayerPara.g_RouteLatQTime.Values[i, j], 1));
                                    pRoute[LagOrder][subID] += dGridOut;
                                    //}
                                }
                                break;
                            case RunOffSourceType.RUNOFF_ELEMENT_BASEQ:
                                {
                                    LagTime = (int)(g_GridLayerPara.g_RouteBaseQTime.Values[i, j] / dTlen);
                                    TimeSpan LagSpan = new TimeSpan(LagTime, 0, 0);
                                    LagOrder = curorder.Add(LagSpan);
                                    //if (LagOrder < totorder)
                                    //{
                                    if (g_GridLayerPara.g_RouteBaseQTime.Values[i, j] <= 0)
                                        g_GridLayerPara.g_RouteBaseQTime.Values[i, j] = 0.1;
                                    dGridOut = snowfactor * snowfactor * GridQ.Values[i, j] / (g_ModelInputPara.BaseQLoss * Math.Pow(g_GridLayerPara.g_RouteBaseQTime.Values[i, j], 1));
                                    pRoute[LagOrder][subID] += dGridOut;
                                    //}
                                }
                                break;
                            default:
                                return false;
                        }
                    }
                }
            }
            return true;
        }

        /// <summary>
        /// 计算植被逐日Albedo
        /// </summary>
        /// <param name="mon"></param>
        /// <param name="day"></param>
        /// <returns></returns>
        public double GetVegAlbedo(int VegID, int mon, int day)
        {
            double result = 0.23;	//取默认值为0.23
            //取真实测量的值
            if (mon >= 1 || mon <= 12)
            {
                result = g_GridLayerPara.VegTypes[VegID].Albedo[mon - 1];
            }
            return result;
        }
        public void LoadSystemParaFile(string path)
        {
            g_ClimatePara.ReadClimateFile(path + "\\SysPara\\Climate.txt");
            g_GridLayerPara.ReadGridIOFile(path);
            //g_GridLayerPara.ReadBasinExtent(m_BasinFile);
            g_ModelRunPara.ReadPara(m_RunParaFile);
            g_ModelInputPara.ReadPara(m_InputParaFile, g_ModelInputPara);
            g_GridResultOut.ReadMidOutFile(m_MidOutFile);
            
        }
        
        /// <summary>
        /// 用于计算河道流速
        /// </summary>
        /// <param name="dUnitQ"></param>
        /// <param name="slp">坡度</param>
        /// <param name="manning_n">曼宁系数</param>
        /// <returns>流速</returns>
        public double ChannelQVelocity(double dUnitQ, double slp, double manning_n)
        {
            double dVflow;
            dVflow = 0.489 * Math.Pow(dUnitQ, 0.25) * Math.Pow(slp, 0.375) / Math.Pow(manning_n, 0.75);
            return dVflow;
        }
        /// <summary>
        /// 计算坡面流速
        /// </summary>
        /// <param name="dUnitQ"></param>
        /// <param name="slp">坡度坡面流速</param>
        /// <param name="manning_n"></param>
        /// <returns></returns>
        public double HillSideQVelocity(double dUnitQ, double slp, double manning_n)
        {
            double dVflow;
            dVflow = Math.Pow(dUnitQ, 0.4) * Math.Pow(slp, 0.3) / Math.Pow(manning_n, 0.6);
            return dVflow;
        }
        /// <summary>
        /// 霍顿下渗曲线计算暴雨过程主函数(刘永和)
        /// </summary>
        public virtual void StormRunoffSim_Horton_Liu()
        {
            int rowCount = g_GridLayerPara.g_DemLayer.rowCount;
            int colCount = g_GridLayerPara.g_DemLayer.colCount;
            this.CellOutSurfFlux = new double[rowCount, colCount];
            this.CellOutLatFlux = new double[rowCount, colCount];
            this.CellOutBaseFlux = new double[rowCount, colCount];
            GridLayerInit();           
            
            //下面这段程序由刘永和所加
            int n = River.CalcuOrderList.Count;
            CurrentSurfQ = new RouteFlux[n];
            CurrentLatQ = new RouteFlux[n];
            CurrentBaseQ = new RouteFlux[n];

            PreSurfQ = new RouteFlux[n];
            PreLatQ = new RouteFlux[n];
            PreBaseQ = new RouteFlux[n];
            for (int i = 0; i < n; i++)
            {
                CurrentSurfQ[i] = new RouteFlux();
                CurrentLatQ[i] = new RouteFlux();
                CurrentBaseQ[i] = new RouteFlux();
                PreSurfQ[i] = new RouteFlux();
                PreLatQ[i] = new RouteFlux();
                PreBaseQ[i] = new RouteFlux();
            }
            //计算间隔时段长度

            StreamWriter sw = new StreamWriter(HydroSimulate.path + "\\Output\\OutletQ.txt");
            sw.WriteLine("Order\tTotalQ\tSurfQ\tLatQ\tBaseQ\tDeepBaseQ");
            //初始化指定出口断面计算结果内存
            this.MyOutletQ = new SortedList<DateTime, OutletQ>();
            //对各降雨时段进行循环
            for (int i = 0; i < g_ClimatePara.pcpdata.saDate.Count; i++)
            {
                DateTime dateTime = g_ClimatePara.pcpdata.saDate[i];
                //取得该降雨时段的dn历
                int dn = DateAndTime.GetDnInYear(dateTime.Year, dateTime.Month, dateTime.Day);
                //对该时段的降雨和气象进行插值
                this.MeteoDataInterp(dateTime);
                //计算该时段的小时数
                double dhr = 1;
                ProgressBar bar = new ProgressBar();
                bar.Text = "正在计算净雨量";
                bar.Show();
                for (int row = 0; row < rowCount; row++)
                {
                    for (int col = 0; col < colCount; col++)
                    {
                        //计算本单元的水量平衡
                        this.CalcCell_Horton(row, col, dhr, dateTime, dn);
                    }
                    bar.progressBar1.Value = (int)(row * 100.0 / rowCount);
                }
                bar.Close();//关闭计算净雨时的bar
                
                //下面开始汇流
                this.CalcRiverRouting_Liu(dateTime, dhr, dn);

                sw.Write(dateTime + "\t");
                sw.Write(this.MyOutletQ[dateTime].m_pOutletQ.ToString("0.0000") + "\t");
                sw.Write(this.MyOutletQ[dateTime].m_pOutletSurfQ.ToString("0.0000") + "\t");
                sw.Write(this.MyOutletQ[dateTime].m_pOutletLatQ.ToString("0.0000") + "\t");
                sw.Write(this.MyOutletQ[dateTime].m_pOutletBaseQ.ToString("0.0000") + "\t");
                sw.Write(this.MyOutletQ[dateTime].m_pOutletDeepBaseQ.ToString("0.0000") + "\t");
                sw.WriteLine();
                sw.Flush();

                MidGridResultOut(dateTime, false);
            }
            sw.Close();
            RiverOutletQ(true);
        }
        /// <summary>
        /// 进行汇流计算
        /// </summary>
        /// <param name="curorder">当前的执行时段序号</param>
        /// <param name="dhr">当前时段的小时数</param>
        /// <param name="totrec">总时段数</param>
        /// <param name="dn">年内的日序</param>
        public void CalcRiverRouting_Liu(DateTime curorder, double dhr, int dn)
        {
            //下面进行坡面漫流汇流
            if (g_ModelRunPara.SurfRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
            {
                RouteFlux[] surfQ = this.MuskCungeSurfRouting_Liu(dhr, g_ModelInputPara.SMCTimeWeight, g_ModelInputPara.SMCGridRTravelTime);

                for (int k = 0; k < River.CalcuOrderList.Count; k++)
                {
                    int row =River.CalcuOrderList[k].row;
                    int col = River.CalcuOrderList[k].col;
                    g_GridResultOut.m_GridRoutingQ[row, col] = surfQ[k].OutFlux;
                }
                //把路径上最后一个像素的出流量作为总出流量
                if (this.MyOutletQ.ContainsKey(curorder) == false)
                    this.MyOutletQ[curorder] = new OutletQ();
                this.MyOutletQ[curorder].m_pOutletSurfQ = surfQ[River.CalcuOrderList.Count - 1].OutFlux;
            }
            else
            {
                this.PureLagGridRouting_Liu(g_GridResultOut.m_GridSurfQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_SURFQ, curorder, 1.0, WaterYearType.WATER_LOW_YEAR);
            }
            //下面进行壤中流汇流
            if (g_ModelRunPara.LatRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
            {
                RouteFlux[] latQ = MuskCungeLatRouting_Liu(dhr,  g_ModelInputPara.LMCTimeWeight, g_ModelInputPara.LMCGridRTravelTime);
                for (int k = 0; k < River.CalcuOrderList.Count; k++)
                {
                    int row = River.CalcuOrderList[k].row;
                    int col = River.CalcuOrderList[k].col;
                    g_GridResultOut.m_GridRoutingQ[row, col] += latQ[k].OutFlux;
                }

                if (this.MyOutletQ.ContainsKey(curorder) == false)
                    this.MyOutletQ[curorder] = new OutletQ();
                this.MyOutletQ[curorder].m_pOutletLatQ = latQ[River.CalcuOrderList.Count - 1].OutFlux;
            }
            else
            {
                PureLagGridRouting_Liu(g_GridResultOut.m_GridLateralQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_LATERALQ, curorder, 1.0, WaterYearType.WATER_LOW_YEAR);
            }
            //下面进行地下径流的汇流
            if (g_ModelRunPara.BaseRouteMethod == SurfRouteMehtod.ROUTE_MUSK_CONGE)
            {
                RouteFlux[] baseQ = this.MuskCungeBaseRouting_Liu(dhr, g_ModelInputPara.BMCTimeWeight, g_ModelInputPara.BMCGridRTravelTime);
                for (int k = 0; k < River.CalcuOrderList.Count; k++)
                {
                    int row = River.CalcuOrderList[k].row;
                    int col = River.CalcuOrderList[k].col;
                    g_GridResultOut.m_GridRoutingQ[row, col] += baseQ[k].OutFlux;
                }

                if (this.MyOutletQ.ContainsKey(curorder) == false)
                    this.MyOutletQ[curorder] = new OutletQ();
                this.MyOutletQ[curorder].m_pOutletBaseQ = baseQ[River.CalcuOrderList.Count - 1].OutFlux;
            }
            else
            {
                PureLagGridRouting_Liu(g_GridResultOut.m_GridBaseQ, dhr, RunOffSourceType.RUNOFF_ELEMENT_BASEQ, curorder, 1.0, WaterYearType.WATER_LOW_YEAR);
            }
            //下面计算深层地下径流的汇流
            if (!this.MyOutletQ.ContainsKey(curorder))
                this.MyOutletQ[curorder] = new OutletQ();
            this.MyOutletQ[curorder].m_pOutletDeepBaseQ = DeepBaseQSim(dn, g_ModelInputPara.DeepBaseQ);
            this.MyOutletQ[curorder].m_pOutletQ = this.MyOutletQ[curorder].m_pOutletSurfQ * g_ModelInputPara.SurfQLinearFactor
                                        + this.MyOutletQ[curorder].m_pOutletLatQ + this.MyOutletQ[curorder].m_pOutletBaseQ + this.MyOutletQ[curorder].m_pOutletDeepBaseQ;
        }
        /// <summary>
        /// 马斯京根-康吉法汇流计算(先合后演)
        /// </summary>
        /// <param name="dT">时间的离散值</param>
        /// <param name="GridQ"></param>
        /// <param name="pRoute"></param>
        /// <param name="pPreRoute"></param>
        /// <param name="x"></param>
        /// <param name="k"></param>
        /// <returns></returns>
        public RouteFlux[] MuskCungeSurfRouting_Liu(double dT, double TimeWeight, double TravelTime)
        {
            ProgressBar bar = new ProgressBar();
            bar.Text = "正在计算汇流";
            bar.Show();
            int n = River.CalcuOrderList.Count;
            for (int i = 0; i < n; i++)
            {
                int row=River.CalcuOrderList[i].row;
                int col=River.CalcuOrderList[i].col;
                //先统计当前栅格(即第i个)的总入流量
                double SumQ = 0;
                for (int index = 0; index < 9; index++)
                {
                    if (index == 4)
                        continue;

                    int R =  row + index / 3 - 1;
                    int C =  col + index % 3 - 1;
                    if (R < 0 || R >= River.demgrid.rowCount)
                        continue;
                    if (C < 0 || C >= River.demgrid.colCount)
                        continue;
                   
                    if (River.m_pDirection[R, C] == 8 - index)
                    {
                        SumQ += this.CellOutSurfFlux[R, C];
                    }
                }
                //if (SumQ < 0)
                //    throw new Exception("");
                CurrentSurfQ[i].InFlux = SumQ;
                //再使用Cunge法计算出流量
                double dX=River.demgrid.resolution;//栅格空间上河流长度值
                if(River.m_pDirection[ row,  col]%2==0)
                    dX*=1.41421;
                
                //如果当前栅格无数据
                if (g_GridResultOut.m_GridSurfQ[row, col] < 0)
                    throw new Exception("流量值为空");
                double tempGridQ = g_GridResultOut.m_GridSurfQ[row, col] * Math.Pow(g_GridResultOut.m_GridSurfQ.resolution, 2) / (1000 * dT * 3600);//栅格面内入流(单位：立方米/秒)
                CurrentSurfQ[i].OutFlux = this.GridRoutingOut_Liu(i, dX, dT, PreSurfQ[i].InFlux, CurrentSurfQ[i].InFlux, PreSurfQ[i].OutFlux, tempGridQ, TimeWeight, TravelTime);
                this.CellOutSurfFlux[row, col] = CurrentSurfQ[i].OutFlux;
                bar.progressBar1.Value = (int)(100.0 * i / n);
            }
            
            //将CurrentQ和PreQ交换
            RouteFlux[] result = CurrentSurfQ;
            RouteFlux[] temp = PreSurfQ;
            PreSurfQ = CurrentSurfQ;
            CurrentSurfQ = temp;
            bar.Close();
            return result;
        }
        public RouteFlux[] MuskCungeLatRouting_Liu(double dT, double TimeWeight, double TravelTime)
        {
            ProgressBar bar = new ProgressBar();
            bar.Text = "正在计算汇流";
            bar.Show();
            int n = River.CalcuOrderList.Count;
            for (int i = 0; i < n; i++)
            {
                int row = River.CalcuOrderList[i].row;
                int col = River.CalcuOrderList[i].col;
                //先统计当前栅格(即第i个)的总入流量
                double SumQ = 0;
                for (int index = 0; index < 9; index++)
                {
                    if (index == 4)
                        continue;

                    int R = row + index / 3 - 1;
                    int C = col + index % 3 - 1;
                    if (R < 0 || R >= River.demgrid.rowCount)
                        continue;
                    if (C < 0 || C >= River.demgrid.colCount)
                        continue;

                    if (River.m_pDirection[R, C] == 8 - index)
                    {
                        if (this.CellOutLatFlux[R, C] < 0)
                            throw new Exception("");
                        SumQ += this.CellOutLatFlux[R, C];
                    }
                }
                if (SumQ < 0)
                    throw new Exception("");
                CurrentLatQ[i].InFlux = SumQ;
                //再使用Cunge法计算出流量
                double dX = River.demgrid.resolution;//栅格空间上河流长度值
                if (River.m_pDirection[row, col] % 2 == 0)
                    dX *= 1.41421;
                
                //如果当前栅格无数据
                if (g_GridResultOut.m_GridLateralQ[row, col] < 0)
                    throw new Exception("流量值为空");
                double tempGridQ = g_GridResultOut.m_GridLateralQ[row, col] * Math.Pow(River.demgrid.resolution, 2) / (1000 * dT * 3600);//栅格面内入流(单位：立方米/秒)
                CurrentLatQ[i].OutFlux = this.GridRoutingOut_Liu(i, dX, dT, PreLatQ[i].InFlux, CurrentLatQ[i].InFlux, PreLatQ[i].OutFlux, tempGridQ, TimeWeight, TravelTime);
                if (CurrentLatQ[i].OutFlux < 0)
                {
                    if (CurrentLatQ[i].OutFlux < -0.1)
                        throw new Exception("");
                    CurrentLatQ[i].OutFlux = 0;
                }
                this.CellOutLatFlux[row, col] = CurrentLatQ[i].OutFlux;
                bar.progressBar1.Value = (int)(100.0 * i / n);
            }

            //将CurrentQ和PreQ交换
            RouteFlux[] result = CurrentLatQ;
            RouteFlux[] temp = PreLatQ;
            PreLatQ = CurrentLatQ;
            CurrentLatQ = temp;
            bar.Close();
            return result;
        }
        public RouteFlux[] MuskCungeBaseRouting_Liu(double dT, double TimeWeight, double TravelTime)
        {
            ProgressBar bar = new ProgressBar();
            bar.Text = "正在计算汇流";
            bar.Show();
            int n = River.CalcuOrderList.Count;
            for (int i = 0; i < n; i++)
            {
                int row = River.CalcuOrderList[i].row;
                int col = River.CalcuOrderList[i].col;
                //先统计当前栅格(即第i个)的总入流量
                double SumQ = 0;
                for (int index = 0; index < 9; index++)
                {
                    if (index == 4)
                        continue;

                    int R = row + index / 3 - 1;
                    int C = col + index % 3 - 1;
                    if (R < 0 || R >= River.demgrid.rowCount)
                        continue;
                    if (C < 0 || C >= River.demgrid.colCount)
                        continue;

                    if (River.m_pDirection[R, C] == 8 - index)
                    {
                        if (this.CellOutBaseFlux[R, C] < 0)
                            throw new Exception("");
                        SumQ += this.CellOutBaseFlux[R, C];
                    }
                }
                if (SumQ < 0)
                    throw new Exception("");
                CurrentBaseQ[i].InFlux = SumQ;
                //再使用Cunge法计算出流量
                double dX = River.demgrid.resolution;//栅格空间上河流长度值
                if (River.m_pDirection[row, col] % 2 == 0)
                    dX *= 1.41421;

                //如果当前栅格无数据
                
                if (g_GridResultOut.m_GridBaseQ[row, col] < 0)
                    throw new Exception("流量值为空");
                double tempGridQ = g_GridResultOut.m_GridBaseQ[row, col] * Math.Pow(River.demgrid.resolution, 2) / (1000 * dT * 3600);//栅格面内入流(单位：立方米/秒)
                CurrentBaseQ[i].OutFlux = this.GridRoutingOut_Liu(i, dX, dT, PreBaseQ[i].InFlux, CurrentBaseQ[i].InFlux, PreBaseQ[i].OutFlux, tempGridQ, TimeWeight, TravelTime);
                if (CurrentBaseQ[i].OutFlux < 0)
                {
                    if (CurrentBaseQ[i].OutFlux < -0.1)
                        throw new Exception("");
                    CurrentBaseQ[i].OutFlux = 0;
                }
                this.CellOutBaseFlux[row, col] = CurrentBaseQ[i].OutFlux;
                bar.progressBar1.Value = (int)(100.0 * i / n);
            }

            //将CurrentQ和PreQ交换
            RouteFlux[] result = CurrentBaseQ;
            RouteFlux[] temp = PreBaseQ;
            PreBaseQ = CurrentBaseQ;
            CurrentBaseQ = temp;
            bar.Close();
            return result;
        }
        /// <summary>
        /// 汇流计算变量
        /// </summary>
        public MuskCungeRoutePara[] RoutePara;
        public MuskingumCunge MuskCunge = new MuskingumCunge();
        public MuskRouteFlux[] pSurfQ;
        public MuskRouteFlux[] pLatQ;
        public MuskRouteFlux[] pBaseQ;
        public double[] pRouteQ;

        public Muskingum Muskingum = new Muskingum();
        public MuskRouteFlux[] pRiverRoute;
        public int m_iNodeNum;
        public double[] m_pX;
        public double[] m_pK;

        /// <summary>
        /// 水量平衡变量
        /// </summary>
        //public GridWaterBalance gridwb;
        //public GreenAmptInfil GreenAmptInfil;
        //public HortonInfil HortonInfil;

        /// <summary>
        /// 土壤水下渗
        /// </summary>
        public static GridLayer m_drateinf;
        public static GridLayer m_drintns;
        public static GridLayer m_dcumr;
        public static GridLayer m_dcuminf;
        public static GridLayer m_dexcum;


        

        //计算结果输出
        public int m_OutRow;
        public int m_OutCol;
        public SortedList<DateTime, OutletQ> MyOutletQ;
        /// <summary>
        /// 与汇流次序表对应的当前流量表(刘永和加)
        /// </summary>
        public RouteFlux[] CurrentSurfQ;
        public RouteFlux[] CurrentLatQ;
        public RouteFlux[] CurrentBaseQ;
        /// <summary>
        /// 与汇流次序表对应的前一时刻流量表(刘永和加)
        /// </summary>
        public RouteFlux[] PreSurfQ;
        public RouteFlux[] PreLatQ;
        public RouteFlux[] PreBaseQ;
        public double[,] CellOutSurfFlux;
        public double[,] CellOutLatFlux;
        public double[,] CellOutBaseFlux;


        //按子流域计算存储
        SortedList<DateTime, double[]> m_pNodeSurfQ;
        SortedList<DateTime, double[]> m_pNodeLatQ;
        SortedList<DateTime, double[]> m_pNodeBaseQ;
        SortedList<DateTime, double[]> m_pNodeOutQ;
        /// <summary>
        /// 这一时段的时长，如24小时或1小时，单位：分钟
        /// </summary>
        public SortedList<DateTime, double> m_pRainInterval;
        int m_iTotalRec;

        //int m_row;
        //int m_col;
        int m_iVegOrd;
        int m_iSoilOrd;


        public bool m_bDate;


        string m_strBeg;
        string m_strEnd;
        /// <summary>
        /// 丰水年或平水年或少水年(其中键为年份值)
        /// </summary>
        SortedList<int, WaterYearType> waterYearType;

        public string m_GridFile;			// GridIO.txt
        public string m_BasinFile;		// Basin.txt 
        public string m_ClimateFile;		// Climate.txt 
        public string m_RunParaFile;		// RunPara.txt 
        public string m_InputParaFile;	// InputPara.txt
        public string m_MuskCoeffFile;	// muskingum_coeff.txt 
        public string m_MidOutFile;		// MidGridOut.txt 
        public string m_WaterYrTypeFile;	// WaterYearType.txt

    }
}